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Abstract 


Interplanetary  shock  waves  (ISWs)  propagating  through  the  solar  wind  can  “collide”  with 
the  earth’s  bow  shock,  resulting  in  a  series  of  new  shocks,  contact  discontinuities,  and 
rarefaction  waves  which  interact  to  effectively  move  the  bow  shock  and  magnetopause 
toward  the  earth.  A  one-dimensional  MacCormack  predictor-corrector  algorithm  with 
Flux  Corrected  Transport  (FCT)  was  developed  to  model  the  ISW-bow  shock  and 
magnetopause  interactions,  and  to  numerically  predict  their  propagation  speeds  after 
collision.  Analytic  relationships  for  the  Mach  numbers  and  propagation  speeds  of  the 
generated  shock  waves  and  contact  discontinuities  were  used  to  validate  the  model 
predicted  propagation  speeds  of  the  moving  bow  shock  to  within  five  percent  of  the 
analytical  solutions.  Propagation  speeds  of  the  moving  magnetopause  were  also 
determined  to  within  five  percent  for  the  gas-dynamic  case. 
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DISPLACEMENT  OF  THE  EARTH’S  BOW  SHOCK  AND  MAGNETOPAUSE  DUE 


TO  AN  IMPINGING  INTERPLANETARY  SHOCK  WAVE 


L  Introduction 


1.1  Problem  Statement 

Predict  and  evaluate  bow  shock  and  magnetopause  motion  resulting  from  the 
interaction  of  an  Interplanetary  Shock  Wave  (ISW)  with  the  earth’s  magnetosheath 
system  using  a  one  dimensional  solution  of  the  gas-dynamic  and  magnetohydrodynamic 
(MHD)  equations. 

1.2  Importance  of  Research 

Rapid,  large  scale  movement  of  the  earth’s  magnetopause  generates  strong 
electric  fields  that  accelerate  charged  particles  to  high  energies,  forming  regions  of 
enhanced  energetic  particle  concentrations  within  the  magnetosphere.  When  the 
displacement  of  the  magnetopause  is  large  enough  to  bring  it  within  approximately  six 
earth  radii  then  satellites  in  geosyncronous  orbits  are  directly  exposed  to  energetic  solar 
wind  particles.  Increased  densities  of  high  energy  particles  in  the  near-earth  environment 
have  damaged  military  and  commercial  satellites,  caused  widespread  radio  and  satellite 
communications  blackouts,  and  have  caused  power  outages  on  earth.  The  ability  to 
predict  magnetopause  motion  resulting  from  the  interaction  of  an  ISW  with  the  earth’s 
magnetosheath  system  would  be  a  valuable  input  for  computer  models  used  to  forecast 
changes  in  ion  and  electron  densities  in  the  near-earth  environment.  Ion  and  electron 
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density  forecasts  based  on  magnetopause  motion  would,  in  turn,  be  used  by  operators  to 
reduce  the  risk  of  damaging  expensive  and  strategically  important  satellite  systems,  and 
by  communicators  to  manage  communications  outages  caused  by  changes  in  the  space 
environment. 

This  thesis  is  especially  relevant  at  this  time  because  of  the  placement  of  the 
Solar  and  Heliographic  Observation  (SOHO)  satellite  into  a  halo  orbit  about  the  LI  point 
upstream  in  the  solar  wind  on  February  14, 1996  (SOHO  home  page,  1997).  From  its 
position  upstream  in  the  solar  wind,  the  SOHO  satellite  is  able  to  provide  density, 
temperature,  and  velocity  measurements  about  an  hour  before  disturbances  in  the  solar 
wind  reach  the  earth.  Using  data  from  SOHO  and  methods  developed  in  this  thesis, 
space  forecasters  will  be  able  to  determine  an  initial  estimate  of  magnetopause  motion 
due  to  ISW  collision  with  the  earth’s  bow  shock  in  advance  of  the  actual  event. 

1.3  Scope  and  Limitations 

This  thesis  will  investigate  the  problem  of  bow  shock  and  magnetopause  motion 
resulting  from  the  interaction  of  an  interplanetary  shock  wave  impinging  on  the  earth’s 
magnetosheath  system.  The  analysis  is  restricted  to  one  dimension.  Magnetic  fields  are 
assumed  to  be  perpendicular  to  the  gas  flow  and  in  the  ^-direction  so  that  slow  and 
intermediate  magnetosonic  shocks  are  eliminated  from  consideration  (Parks,  1991,418). 
The  difficulties  associated  with  a  kinetic  approach  are  avoided  in  favor  of  two  fluid 
models:  a  gas-dynamic,  ideal  fluid  approximation  in  which  the  magnetic  field  is 
neglected,  and  a  magnetohydrodynamic  fluid  approximation  where  the  effects  of 
magnetic  and  electric  fields  on  the  flow  are  taken  into  account.  A  one-dimensional 
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treatment  can  not  seamlessly  simulate  the  entire  interaction  process  from  ISW  collision 
with  the  bow  shock  to  magnetopause  displacement;  consequently,  this  thesis  will  explore 
the  ISW-bow  shock  interaction  and  magnetopause  motion  as  two  separate  steps. 

1.4  Research  Objectives 

This  thesis  has  four  objectives.  The  first  objective  is  to  develop  an  understanding 
of  the  physical  processes  associated  with  the  displacement  of  the  earth’s  bow  shock  and 
magnetopause  initiated  by  collision  of  an  ISW  with  the  magnetosheath  system.  The 
second  objective  is  to  develop  a  one-dimensional  numerical  model  of  the  ISW- 
magnetosheath  interactions  based  upon  gas-dynamic  and  MHD  fluid  equations  that  will 
allow  determination  of  the  displacements  and  velocities  of  the  bow  shock  and  the 
magnetopause.  A  third  objective  is  to  review  and  understand  prior  analytic  treatments  of 
the  ISW-magnetosheath  interactions,  providing  a  basis  for  comparison  with  numerical 
results.  And  the  fourth  objective  is  to  evaluate  the  performance  of  the  numerical  model 
by  comparing  simulation  results  with  analytic  solutions. 

1.5  Thesis  Overview 

Following  this  introduction.  Chapter  II  will  provide  a  general  background  on 
plasma  and  shock  physics.  The  background  chapter  begins  by  citing  the  March  24,  1991 
Storm  Sudden  Commencement  (SSC)  event  that  resulted  from  an  ISW  impinging  upon 
the  earth’s  bow  shock.  A  discussion  of  the  near-earth  environment  is  given  after  the  SSC 
reference  to  establish  the  physical  setting  of  the  problem.  Then,  in  the  next  section,  some 
basic  ideas  about  shock  waves  in  ordinary  gases  and  in  plasmas  are  presented,  along  with 
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a  discussion  of  the  Rankine-Hugoniot  relations  that  describe  the  conservation  of  mass 
momentum,  and  energy  across  a  shock  boundary.  The  discussion  of  shock  waves  is 
followed  by  a  brief  examination  of  the  interactions  of  shock  waves  with  other  shock 
waves,  and  shock  waves  with  regions  of  discontinuous  jumps  in  density  and  temperature 
within  a  gas  or  plasma.  Chapter  n  ends  with  a  summary  of  the  background  material 
covered  to  that  point. 

Chapter  III,  a  separate  literature  review,  was  deemed  necessary  to  discuss  in  detail 
three  articles  that  specifically  examine  the  interaction  of  an  interplanetary  shock  wave 
with  the  earth’s  magnetosheath  system.  The  first  article,  by  W.  W.  Shen  and  M.  Dryer 
(1972),  considers  the  interaction  of  an  interplanetary  shock  wave  with  the  bow  shock;  the 
second  article,  by  Shen  (1973)  alone,  describes  the  interaction  of  a  shock  wave  with  the 
magnetopause.  The  third  article  included  in  the  literature  review  is  a  detailed  work  by 
S.  A.  Gnb,  et.  al.  (1979)  which  develops  analytic  solutions  for  shock  wave-bow  shock 
and  shock  wave-magnetopause  interactions  in  both  the  gas-dynamic  and  MHD  cases. 

The  solutions  developed  by  Grib  provide  a  means  to  analytically  determine  position  and 
velocity  for  the  displaced  bow  shock  and  magnetopause  as  a  function  of  the  Mach 
number  of  the  incident  interplanetary  shock  wave;  they  also  serve  as  a  comparison  to  test 
how  well  the  numerical  simulation  predicts  bow  shock  and  magnetopause  motion. 

Methodology  is  the  subject  of  Chapter  IV.  The  first  section  is  a  presentation  of 
the  design  and  development  of  a  numerical  model  to  simulate  the  ISW  interaction  with 
the  earth’s  magnetosheath  system.  Development  of  the  numerical  model  begins  vvith  the 
fluid  conservation  equations,  from  which  a  set  of  time  and  spatial  dependent  differential 
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equations  for  the  conservation  variables — i.e.  mass  density,  momentum,  magnetic  field, 
and  total  energy — are  derived.  Next,  the  numerical  algorithm  for  a  one-dimensional, 
explicit,  time  dependent  solution  is  discussed.  This  section  also  presents  the  procedure 
used  to  extract  values  for  pressure,  density,  temperature,  magnetic  field,  and  flow 
velocity  from  updated  conservation  variables,  as  well  as  the  method  for  determining 
shock  front  and  contact  discontinuity  locations  at  each  time  step.  After  detailing  the 
development  of  a  numerical  model,  the  next  section  examines  the  methodology  of  model 
validation. 

Chapter  V  is  a  presentation  of  the  simulation  for  four  different  ISW  Mach 
numbers  in  both  the  gas-dynamic  and  the  MHD  cases.  Here  numerical  solutions  for  bow 
shock  and  magnetopause  position  and  velocity  are  given.  Bow  shock  and  magnetopause 
position  and  velocity  determined  from  solutions  to  the  analytic  equations  are  also 
presented  in  this  chapter,  and  the  two  sets  of  solutions  are  quantitatively  compared. 

Based  on  the  results  from  Chapter  V,  the  final  chapter  contains  conclusions  about 
model  performance,  its  strengths  and  weaknesses,  and  its  applicability  to  the  problem  of 
determining  bow  shock  and  magnetopause  displacement  due  to  an  impinging  ISW. 

Some  comments  are  made  concerning  accomplishment  of  the  research  objectives  and 
extension  of  this  work  beyond  an  investigation  into  separate  bow  shock  and 
magnetopause  movements  to  a  continuous  treatment  of  the  entire  magnetosheath 
displacement.  And  finally,  at  the  end  of  this  chapter,  recommendations  regarding  model 
improvements  and  real-world  validation  are  presented. 
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n.  Background 


At  03:41  UT  on  March  24, 1991  the  formation  of  a  new  radiation  belt  near  the 
equatorial  plane  at  2.55  earth  radii  was  observed  by  the  Combined  Release  and  Radiation 
Effect  Satellite  (CRRES)  at  a  point  in  its  orbit  which  coincidentally  passed  through  the 
inner  edge  of  the  formation  region.  This  new  radiation  belt  resulted  from  the  interaction 
of  an  interplanetary  shock  wave  which  compressed  the  earth’s  magnetopause  inside 
geostationary  orbit  and  caused  a  surge  in  the  geomagnetic  field  characteristic  of  a  Storm 
Sudden  Conunencement  (SSC)  event  (Li„  1993:2423).  Magnetopause  compression  was 
accompanied  by  electron  and  ion  drift  echo  events  from  which  it  was  determined  that 
energies  of  the  electrons  injected  into  this  new  radiation  belt  were  greater  than  15  MeV 
(Li,  1993:2423),  and  that  the  injected  proton  energies  were  in  the  range  20-80  MeV 
(Hudson,  1995:291).  The  radiation  belt  formed  in  less  than  150  seconds  and  persisted 
beyond  the  end  of  the  CRRES  mission  six  months  later  (Li,  1993:2423).  Xinlin  Li 
successfully  modeled  the  formation  of  the  electron  component  of  the  new  radiation  belt 
resulting  from  the  March  1991  SCC  using  a  relativistic  guiding  center  code  (Li, 
1993:2423);  later  M.K.  Hudson  used  the  same  technique  to  model  the  proton  component 
(Hudson,  1995:291). 

The  preceding  example  is  cited  to  illustrate  three  ideas:  that  interplanetary  shock 
waves,  when  they  interact  with  the  earth’s  magnetosphere,  are  capable  of  producing 
major  changes  in  the  near-earth  enviromnent;  that  compression,  or  motion,  of  the  earth’s 
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magnetopause  is  associated  with  changes  in  particle  energies  and  distributions  within  the 
magnetosphere;  and  that  it  is  possible  to  simulate  these  changes  using  a  numerical 
computer  model. 

2.1  Space  Environment 

To  understanding  how  interactions  with  interplanetary  shock  waves  produce 
motions  of  the  bow  shock  and  magnetopause  it  is  first  necessary  to  understand  the 
physical  characteristics  of  the  near-earth  environment.  Figure  1  is  a  representation  of  the 
near-earth  enviromnent. 


Bow  Shock 


Figure  1.  Earth  bow  shock-magnetopause  system.  Scale  lengths  are  taken 
from  Parks  (1996:501-502). 
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Although  it  is  commonly  assumed  that  space  is  a  vacuum,  within  the  solar  system 
there  is  in  fact  a  tenuous  gas  made  up  almost  entirely  of  ionized  hydrogen  atoms  in  the 
form  of  free  protons  and  electrons  streaming  outward  from  the  sun  at  velocities  greater 
than  the  local  speed  of  soimd  within  the  gas.  This  streaming  ionized  gas,  or  plasma, 
called  the  solar  wind  not  only  interacts  with  other  matter  kinetically  as  ordinary  gases  do, 
but  its  motion  is  also  influenced  by  electric  and  magnetic  fields  in  space.  The  properties 
of  the  solar  wind  vary,  with  values  of  flow  speed,  density,  and  magnetic  field  linked  to 
changes  on  the  sun’s  surface.  Table  1  lists  typical  solar  wind  parameters,  taken  from 
Grib  (1979:5909),  that  are  used  in  numeric  and  analytic  calculations  later  in  this  thesis. 


Table  1.  Solar  Wind  Parameters 


Parameter 

Value 

Density ,  po  (kg/m^) 

1.837x10'^" 

Temperature,  To  (®K) 

4.000  xlO'^ 

Flow  Velocity,  Uo  (km/sec) 

280.0 

Magnetic  field.  Bo  (Tesla) 

3.500  xlO’^ 

Pressure,  Po  (Pa) 

0.121  xlO''" 

Sound  Speed,  ao  (km/sec) 

33.20 

Thermal  Velocity,  Vth  (km/sec) 

18.18 

Plasma  Beta,  Po 

2.490 

As  the  solar  wind  flows  toward  the  earth  it  encounters  the  magnetosphere,  a 
region  where  the  relatively  strong  magnetic  field  of  the  earth  forms  a  protective  cavity 
around  which  the  solar  wind  is  constrained  to  flow  imder  the  influence  of 
electromagnetic  forces.  Upon  encountering  the  magnetosphere,  a  bow  shock  wave  forms 
in  front  of  the  magnetopause  similar  to  the  aerodynamic  shock  found  in  front  of  a  blunt 
obstacle  in  supersonic  wind  tuimel  experiments  (Tascione,  1994:57).  At  the  bow  shock. 
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solar  wind  plasma  undergoes  rapid  and  dramatic  increases  in  pressure,  density, 
temperature,  and  magnetic  field  across  the  span  of  a  few  hundred  kilometers;  flow  speed 
across  the  bow  shock  also  changes  from  supersonic  to  subsonic  flow  (Grib:  1979:5909). 
From  data  presented  by  Parks  (1996:501),  the  bow  shock  appears  to  maintain  a  fairly 
consistent  stand-off  distance  from  the  magnetopause  of  about  two  earth  radii. 

Between  the  bow  shock  and  the  magnetopause  lies  a  region  of  compressed  and 
heated  solar  wind  plasma  called  the  magnetosheath  (Tascione,  1994:61).  Here  the 
magnetic  field  is  more  disordered  than  in  the  solar  wind  or  in  the  magnetosphere,  and  the 
plasma  is  irregularly  distributed  through  the  region  (Tascione,  1994:57). 

The  magnetopause  is  the  boundary  separating  the  interplanetary  medium  from  the 
magnetosphere  and  is  the  surface  where  the  outward  force  of  the  total  magnetospheric 
pressure  is  balance  by  the  force  of  the  solar  wind  plasma  pressure  (Tascione,  1994:57). 
Total  pressure,  p*,  is  defined  by  Eq  (1) 


where, 

p  =  gas  dynamic  pressure, 

and  the  second  term  on  the  right  is  magnetic  pressure  with, 
B  =  magnitude  of  the  magnetic  field, 

Po  =  permeability  of  free  space. 


(1) 
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The  magnetopause  is  classified  as  a  tangential  discontinuity  (Parks,1996;335) 
meaning  that,  at  the  boundary,  plasma  flow  velocity  and  magnetic  field  components 
normal  to  the  surface  are  zero — ^they  do  not  penetrate  into  the  magnetosphere;  density, 
temperature,  and  magnetic  field,  however,  are  allowed  to  change  across  the  boundary  to 
preserve  the  continuity  of  total  pressures  (Parks,  1996:329).  Classifying  the 
magnetopause  as  a  tangential  discontinuity  requires  that  the  solar  wind  decelerate 
through  the  magnetosheath  so  that,  at  the  magnetopause,  the  velocity  component  normal 
to  the  surface  satisfies  the  boxmdary  condition — ^the  point  at  which  the  velocity 
component  becomes  zero  is  called  the  stagnation  point  (Parks,  1996:317).  The  position 
of  the  magnetopause  is  very  sensitive  to  changes  in  solar  wind  density,  flow  velocity,  and 
to  changes  in  strength  and  orientation  of  the  Interplanetary  Magnetic  Field  (IMF) 
(Tascione,  1996:57). 

Inside  the  magnetopause  lies  the  magnetosphere — ^the  region  of  space 
immediately  surrounding  the  earth  permeated  solely  by  the  earth’s  magnetic  field.  It  can 
be  shown  that  the  magnetic  flux  through  a  volume  of  solar  wind  plasma  remains 
essentially  unchanged,  or  “frozen-in”,  so  that  the  IMF  carried  along  by  the  solar  wind 
does  not  penetrate  into  the  magnetosphere  as  the  flow  is  forced  around  the 
magnetopause.  Satellite  observations  have  shown  that  plasma  density  is  lower,  and 
temperature  higher,  just  inside  the  magnetosphere  than  in  the  magnetosheath  (Parks, 
1996:338-39). 
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2.2  Shocks 


A  discussion  of  the  earth’s  bow  shock  leads  naturally  to  a  closer  examination  of 
shock  waves  and  shocks  in  space. 

2.2.1  Sound  Speed  in  Gases.  In  ordinary  gases,  collisions  between  gas  particles 
transfer  momentum  and  energy  among  the  molecules,  allowing  compression  acoustic 
waves  to  propagate  through  the  medium.  In  an  ideal  gas,  neglecting  viscosity  and  heat 
conduction,  and  assuming  an  adiabatic  compression,  a  disturbance  will  propagate 
through  the  gas  at  the  speed  of  sound.  Thus  the  speed  of  sound  can  be  thought  of  as  the 
speed  at  which  information  about  a  change  in  the  state  of  a  gas,  or  perturbation,  is 
transmitted  (Kivelson  &  Russell,1995;  130-131)  and  is  given  for  an  ideal  isentropic  gas 
by. 


where, 

a  =  speed  of  sound, 

Y  =  isentropic  exponent  s  ratio  of  specific  heats  Cp  /  Cy,  =  5/3  for  an  ideal, 

V  monatomic  gas, 
p  =  gas  pressure, 
p  =  mass  density. 

2.2.2  Shock  Formation.  Shocks  are  formed  by  the  steepening  of  compression 
waves  in  a  gas.  The  speed  of  sound  is  the  speed  limit  of  an  ordinary  wave  in  a  gas:  a 
disturbance  can  not  propagate  through  a  gas  faster  than  the  sound  speed  unless  the  gas 
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changes  from  its  original  state  in  such  a  way  that  the  sound  speed  increases.  If  the  gas  is 
described  by  the  isentropic  form  of  the  equation  of  state,  p/p^  =  constant  (Wright, 
1961:6);  then  for  a  disturbance  resulting  in  a  compression  wave  it  can  be  shown  that  the 
new  speed  of  sound,  ai,  is 


(3) 

where, 

ao  =  sound  speed  given  by  Eq  (2)  for  po,  po  in  the  undisturbed  flow. 

Pi  =  higher  density  within  the  compression  wave. 

Density  pi  is  greater  within  the  compression  wave  than  in  the  undisturbed  gas;  ai  is 
therefore  greater  than  ao-  Because  the  speed  of  sound  is  greatest  at  the  peak  of  the 
compression  wave  where  the  density  is  highest,  the  peak  of  the  wave  will  propagate 
faster  and  catch  up  to  the  front  of  the  wave  ahead  of  it,  steepening  the  perturbation 
gradient  of  the  wave  imtil  the  flow  becomes  non-adiabatic,  i.e.  energy  begins  to  flow  into 
or  out  of  the  system  so  that  changes  become  irreversible.  At  that  point  a  balance  is 
reached  between  steepening  and  dissipation  due  to  viscosity  and  thermal  conduction, 
resulting  in  a  stable  form  of  the  wave  front  called  a  shock  (Parks,  1996:416). 

To  illustrate  the  balance  between  steepening  and  dissipation,  and  their  role  in 
shock  formation,  an  analytic  solution  to  Burger's  Equation  is  discussed.  Burger’s 
Equation  describes  the  propagation  of  a  perturbation  wave  in  a  non-linear  fluid  and  is 
given  by. 
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(4) 


where, 

p(x,  t)  =  a  perturbation  in  a  fluid  as  a  flmction  of  position,  x,  and  time,  t; 
c  =  propagation  speed  of  the  wave; 


a  =  coefficient  of  the  non-linear  term  p(x,t) 


^P(x,t) 

dx 


a  =  wave  dampening  coefficient  (Landau,  1959:351). 

Transforming  to  a  wave-centered  reference  frame  moving  along  with  the  wave  at  speed  c, 
Eq  (  4 )  becomes. 


here. 


dp'(^,t)  ..ap'(^,t)  5p'(^,t) 


(5) 


p'(4,t)  =  -  ap(x,t) 

4  =  characteristic  of  the  wave  propagation,  in  this  case  a  wave  moving  to  the 
right  given  by  x  +  c  t, 
p  =  dissipation  term,  a  c^. 

The  solution  to  the  transformed  Burger’s  Equation  (Eq  (5)),  p'  (^,  t),  is  derived  using  the 
techniques  of  Vvedensky  (1992:274-79,361-66)  and  has  the  form. 
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From  Eq  ( 6 )  steepening  and  broadening  of  the  wave  front  depend  on  the  values 
of  a  and  p:  a  larger  coefficient  of  non-linearity  causes  the  wave  to  steepen  faster  in  time 
and  to  have  a  steeper  slope;  greater  dissipation  causes  the  wave  to  steepen  more  slowly 
and  to  have  a  more  inclined  slope.  When  non-linear  and  dissipative  terms  in  Eq  ( 5  )  are 
equal,  then  5p(^,  t)/5t  =  0  and  a  steady  wave  form  is  achieved. 

The  same  effects  of  wave  steepening  manifest  by  Eq  (  6  )  are  observed  when  a 
disturbance  propagates  through  a  gas  at  supersonic  speeds.  At  speeds  greater  than  the 
sound  speed  non-linear  effects  proportional  to  the  velocity  become  important  so  that  a 
wave  front  steepens,  and  when  dissipation  within  the  gas  balances  steepening,  a  steady 
shock  wave  is  formed  For  pinposes  of  this  thesis  a  shock  is  defined  as  a  disturbance 
propagating  at  supersonic  speed  through  a  plasma  as  a  steady  wave  front  that  irreversibly 
alters  the  state  of  the  gas  behind  it.  (Kivelson  &  Russell,  1995: 131) 

2.2.3  Shocks  in  Space.  Shocks  in  space,  however,  are  more  complicated 
because  interplanetary  plasmas  are  more  rarefied  than  ordinary  gases  and  the  presence  of 
electric  and  magnetic  fields  in  space  affect  plasma  wave  propagation.  According  to 
Kivelson  and  Russell  (1995:129),  “Collisions  in  an  ordinary  gas  serve  to  transfer 
momentum  and  energy  among  the  molecules,  and  they  provide  the  coupling  that  allows 
the  basic  wave,  the  sound  wave,  to  exist.”  Plasmas  in  space  are  so  tenuous  that 
collisional  coupling  is  essentially  absent  (Kivelson  and  Russell,  1995:129).  Parks 
(1996:414)  indicates  that  the  exact  mechanisms  for  shock  formation  in  space  plasmas  are 
not  yet  identified  but  suggests  that  plasma  collective  effects  and  wave-particle 
interactions  may  play  a  part  in  the  processes.  If  in  the  macroscopic  regime  the  following 
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assumptions  are  made:  that  the  charged  plasma  particles  behave  as  a  single  entity  and 
that  Coulomb  forces  between  them  are  negligible,  the  plasma  is  in  thermal  equilibrium, 
and  the  motion  of  the  plasma  is  affected  by  electric  and  magnetic  fields;  then  the  plasma 
behaves  as  an  ideal  fluid~this  is  the  Magnetohydrodynamic  (MHD)  approximation 
(Parks,  1996:141).  Plasmas  behaving  as  ideal  fluids  in  the  MHD  approximation  are 
governed  by  the  same  relationships  that  allow  shock  waves  to  form  in  ordinary  gases;  the 
difference  is  that  in  the  MHD  case  there  are  four  different  types  of  waves,  three  of  which 
are  compressional  and  can  lead  to  shock  formation  in  space. 

2«2.4.  Magnetosonic  Shocks.  The  four  types  of  waves  possible  in  a  magnetized 
plasma  are  the  acoustic  wave,  the  Alfven  wave,  and  the  fast  and  slow  magnetosonic 
waves — of  these  four,  only  the  Alfven  wave  is  not  compressional.  For  propagation 
perpendicular  to  the  magnetic  field,  the  AlfVen  and  slow  magnetosonic  waves  vanish 
(Parks,  1996:41 8)  so  that  only  acoustic  and  fast  magnetosonic  waves  are  available  to 
form  shocks.  Just  like  the  acoustic  wave,  the  fast  magnetosonic  wave  has  a  propagation 
speed,  called  the  fast  magnetosonic  sound  speed,  which  depends  on  the  characteristics  of 
the  plasma  through  which  it  travels.  The  fast  magnetosonic  sound  speed,  Cf ,  is  given  by 

c,=Va'+V/  (7) 

where, 

a  =  sound  speed  (Eq  ( 1 )), 

Va  s  Alfven  wave  speed  =  (B^po  p)*^, 

B  =  magnetic  field. 
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Similar  to  an  acoustic  wave,  if  a  magnetosonic  wave  propagates  through  a  magnetized 
plasma  at  a  speed  greater  than  Cf ,  then  it,  too,  will  steepen  to  become  a  fast 
magnetosonic  shock.  The  fast  magnetosonic  speed  of  magnetospheric  plasma  is  a  key 
parameter  in  determining  magnetopause  motion. 

2,2.5  Rankine-Hueoniot  Relations.  As  discussed  in  Section  2.2.2,  shock  waves 
are  non-adiabatic  regions  where  kinetic  energy  of  the  fluid  is  converted  to  thermal  energy 
through  dissipation;  the  total  energy  of  the  fluid,  however,  must  remain  the  same. 
Conversion  between  forms  of  energy  at  a  shock  necessitates  changes  in  plasma 
parameters.  Figure  3  is  representative  of  changes,  or  jumps,  in  plasma  parameters  across 
a  generic  shock  located  at  the  mid-point  of  each  panel;  numbers  in  each  panel  indicate 
unshocked,  0,  and  shocked,  1,  plasma. 
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Figure  3.  Shock  jump  conditions 
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Conservation  of  mass,  momentum,  and  energy  across  a  shock  require  the 


following  relations  to  be  true; 


[pu]  =  0 

(8) 

[pu^  +  p]  =  0 

(9) 

2  ^  ^  1 

pu  YP  1 

=  0 

(10) 

2  (Y-l)J" 

where,  u  =  gas  flow  velocity  (Parks,  1996:421).  The  bracket  notation  indicates  the 
difference  between  the  included  quantities  on  the  two  sides  of  a  shock;  e.g.,  [p  u  ]  =  0 
means  pi  Ui  -  po  Uo  =  0,  where  the  subscripts  1,  and  0  correspond  to  the  shocked,  and 
unshocked  regions  of  Figure  3.  Equations  (  8  )  -  ( 10 )  are  the  Rankine-Hugoniot 
relations  in  the  gas-dynamic  case  and  from  them  the  ratios  of  density,  pressure,  flow 
velocity,  and  temperature  between  the  two  sides  of  a  shock  are  determined.  These  ratios, 
or  jump  conditions,  are  valid  in  the  shock-centered  reference  frame  and  are  given  by. 


Pi_^ 

(Y  +  1)M^ 

(11) 

Po 

(Y-1)M^+2  “  ^ 

Pi__ 

2yM^-(Y-1) 

(12) 

Po 

(r+i) 

iiL_ 

(T-1)M"+2 

(13) 

Uo 

(Y-i-I)M^  ^ 

Il_  y 

To  T| 

(14) 
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where, 

M  s  Mach  number  =  ratio  of  the  shock  propagation  speed,  V,  relative  to  the 
upstream  speed  of  soimd  =  V/ao, 

T|  =  shock  compression  ratio,  or  compressibility, 
y  =  shock  strength  (Parks,  1991:422). 


From  Equations  (15)  -  (18)  the  MHD  jump  conditions  are  found: 


(19) 

(20) 

(21) 

(22) 
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^  s  degree  of  compressibility  =  ti*‘  (Eq(13)), 

Po  “  plasma  beta  of  the  upstream  gas  =  2  |4o  Pc/Bo^, 

p  =  (y-1)/(y+1)  (Grib,  1979:5908).  Here  again,  subscripts  refer  to  the  shocked  and 
imshocked  regions  of  Figure  3.  When  the  magnetic  field  is  zero  P  goes  to  infinity  and  the 
MHD  jump  conditions  reduce  to  the  gas-dynamic  jump  conditions.  From  Eq  (20)  the 
magnetic  field  ratio  is  equal  to  C  meaning  that  B  varies  as  density  across  a  shock. 

Using  the  jump  conditions,  plasma  parameters  on  either  side  of  a  shock  can  be  found  if 
the  shock  Mach  number  and  the  parameters  on  the  other  side  are  known. 

2.3  Shock  Interactions 

After  discussing  the  near-earth  environment  and  shocks,  the  final  piece  of 
background  information  presented  is  an  examination  of  interactions  between  shocks, 
discontinuities,  and  rafefaction  waves — waves  in  which  the  density  and  pressure 
decrease,  while  the  magnitude  of  the  velocity  increases.  An  incoming  interplanetary 
shock  wave  (ISW)  is  a  shock  wave  propagating  through  the  ambient  solar  toward  the 
earth.  As  it  approaches  earth,  the  ISW  interacts  with  the  bow  shock — a  shock-shock 
interaction — generating  two  new  shock  waves  that  move  inward  toward  the  earth  at 
different  speeds.  The  slower  shock  wave  in  effect  becomes  the  new  bow  shock;  the 
faster  shock  wave  rapidly  traverses  the  magnetosheath  to  interact  with  the 
magnetopause — a  shock-discontinuity  interaction.  The  interaction  of  the  fast  shock  wave 


with  the  magnetopause  creates  a  third  new  shock  wave  transmitted  into  the 
magnetosphere  and  a  new  contact  discontinuity  that  becomes  a  new  magnetopause  also 
moving  toward  the  earth.  Shock-shock  and  shock-discontinuity  interactions  will  be 
examined  separately  in  the  next  two  sections. 

2.3.1  Shock-Shock  Interactions.  According  to  Landau  and  Lifshitz,  when  two 
shock  waves  intersect  the  result  will  be  two  new  shock  waves  separated  by  a  contact 
discontinuity  (1959:408-409).  A  contact  discontinuity  is  a  surface  characterized  by 
jumps  in  density  and  temperature  from  one  side  to  the  other,  but  across  which  plasma 
velocity  and  pressure  are  continuous.  The  following  continuity  relationships  therefore 
apply  across  a  contact  discontinuity: 


[u]  =  0 

(24) 

[p]^0 

(25) 

p+- —  =0 

2  jIq 

(26) 

Courant  and  Friedrichs  (1948: 179)  reached  a  similar  conclusion  when  they  determined 
that  collisions  between  two  shocks  of  different  strengths  produces  two  new  shock  waves 
moving  away  from  each  other  in  a  reference  frame  moving  with  the  contact  discontinuity. 

Figure  4  depicts  a  shock-shock  interaction  representing  the  ISW-bow  shock 
collision  on  an  x-t  (distance-time)  diagram.  In  this  figure  the  ISW  and  bow  shock 
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intersect  at  time  tj  resulting  in  two  new  shocks,  S3  and  S4,  with  a  contact  discontinuity. 
Cl,  between  them  consistent  with  Landau  and  Lifshitz. 


X  (distance) 

Figure  4.  x-t  diagram  for  shock-shock  interaction 


The  propagation  speed  of  a  shock  is  determined  from  an  x-t  diagram  by  the  inverse  of  the 
slope  drawn  through  its  x,  t  coordinates;  faster  shocks  have  a  slope  that  is  more 
horizontal  while  a  vertical  line  represents  a  stationary  feature.  The  reference  to  Courant 
and  Friedrichs  would  seem  to  indicate  that  Ci  in  Figure  4  should  be  drawn  as  a  vertical 
line  with  S3  and  S4  propagating  in  opposite  directions;  that  all  three  structures  are  shown 
traveling  to  the  right  in  the  positive  x  direction  is  an  illustration  that  they  propagate  along 
with  the  bulk  plasma  flow  toward  the  earth. 


22 


The  numbers  in  Figure  4  correspond  to  regions  of  plasma  “shocked”,  or  altered  by 
one  or  more  of  the  depicted  shock  waves.  By  convention,  in  a  reference  frame  traveling 
with  a  shock  wave,  unshocked  plasma  flowing  into  the  shock  is  designated  as  the 
“upstream”  plasma,  while  “shocked”  plasma  flowing  out  is  referred  to  as  the 
“downstream”  plasma.  Region  0  refers  to  the  pre-interaction,  unshocked,  ambient  solar 
wind,  while  Region  1  represents  plasma  from  Region  0  shocked  by  the  stationary  bow 
shock.  In  a  similar  manner  Region  2  depicts  upstream  solar  wind  plasma  shocked  by  an 
advancing  ISW.  After  ti,  the  upstream  plasma  shown  flowing  into  S3  comes  from  Region 
2,  and  Region  1  plasma  becomes  the  upstream  flow  into  S4;  3  and  4  therefore  depict 
regions  of  twice-shocked  plasma. 

2.3.2  Shock-Discontinuity  Interactions.  Turning  again  to  Landau  and  Lifshitz, 
(1959:408)  the  collision  between  a  shock  wave  and  a  tangential  discontinuity  produces  a 
reflected  rarefaction  wave,  a  new  shock  wave,  and  a  contact  discontinuity.  Courant  and 
Friedrichs  (1948: 179)  establish  that  if  the  second  medium  has  a  lower  sound  speed,  or  if 
the  incident  shock  is  sufficiently  weak,  then  a  rarefaction-shock  wave  pair  is  produced 
when  a  shock  impinges  on  a  contact  discontinuity.  Figure  5  is  an  x-t  depiction  of  the 
interaction  between  the  shock  wave  generated  by  the  ISW-bow  shock  intersection,  S4, 
and  the  magnetopause  represented  as  a  tangential  discontinuity. 


23 


A 


X  (distance) 


Figure  5.  S4  -  magnetopause  interaction 

Here,  the  time  is  extended  so  that  S4  collides  with  the  magnetopause  at  ta  to  produce  a 
shock  transmitted  into  the  magnetosphere,S6,  a  reflected  rarefaction  wave,  R5,  and  a 
contact  discontinuity  between  them.  Cm,  representing  the  new  magnetopause  moving 
toward  the  earth.  Regions  1, 2, 3,  and  4  are  the  same  as  in  Figure  4.  Regions  5  and  6  are 
the  regions  downstream  from  R5,  Se,  and  Region  m  is  the  undisturbed  magnetosphere. 

The  slope  of  Se  was  drawn  at  a  more  horizontal  angle  to  the  x-axis  to  illustrate 
that,  upon  crossing  into  a  region  of  lower  density  and  higher  temperature,  the  propagation 
speed  of  a  shock  wave  increases  Wright  (1961 :79).  A  second  point  regarding  Figure  5  is 
that,  although  R5  is  depicted  propagating  to  the  left,  this  may  not  be  necessarily  true.  A 
rarefaction  wave  is  not  a  shock  and  thus  propagates  at  the  local  speed  of  sound  within  the 
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plasma.  If  the  soimd  speed  is  less  than  the  flow  speed  of  the  plasma  in  Region  4,  then  R5 
will  propagate  along  with  the  plasma  flow  to  the  right  (Grib,  1979:591 1). 

2.4  Background  Summary 

Before  moving  on  to  the  next  chapter  a  brief  summary  of  the  background  material 
may  be  helpful.  The  author  can  only  handle  simple  ideas,  so  this  summary  will  attempt 
to  distill  the  important  points  into  nut-shell  concepts  that  the  reader  can  take  with  him 
through  the  remainder  of  this  thesis. 

Three  broad  topics  were  touched  upon:  the  near-earth  environment,  shocks  in 
plasmas,  and  interactions  between  shocks.  The  near-earth  environment  is  “filled”  with  a 
very  thin  plasma  flowing  supersonically  toward  the  earth  from  the  sun.  The  formation  of 
shock  waves  in  space  plasmas  are  more  difficult  to  understand,  but  if  a  few  simplifying 
assumptions  about  how  plasma  particles  interact  with  each  other  and  with  electomagnetic 
fields,  then  to  a  good  approximation  the  plasma  behaves  as  an  ideal  gas.  Shocks  are 
boimdaries  across  which  plasma  parameters  change  rapidly.  Fortunately  a  set  of  simple 
mathematical  relationships  exist  which  allow  plasma  parameters  on  one  side  of  a  shock 
or  discontinuity  to  be  determined  if  parameters  on  the  other  side  are  known.  And  finally, 
in  the  case  of  an  interplanetary  shock  wave  impinging  on  the  magnetosheath  system,  bow 
shock  and  magnetopause  displacement  is  determined  by  shock-shock  and  shock- 
discontinuity  interactions.  The  sequence  of  events  is  as  follows:  first,  an  ISW  impinges 
upon  the  stationary  bow  shock  creating  two  new  shock  waves:  a  fast  shock  which  quickly 
traverses  the  magnetosheath  to  interact  with  the  magnetopause,  and  a  slower  shock  that 
becomes  the  new  bow  shock  moving  in  towards  the  earth;  and  then,  at  the  magnetopause. 
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the  fast  shock  wave  collides,  resulting  is  a  third  shock  wave  transmitted  into  the 
magnetosphere,  a  rarefaction  wave  propagating  back  toward  the  advancing  bow  shock  (if 
the  plasma  flow  speed  in  the  disturbed  magnetosheath  is  less  than  the  speed  of  sound), 
and  a  new  magnetopause  also  moving  in  toward  the  earth. 
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in.  Literature  Review 


A  review  of  the  pertinent  literature  yielded  three  articles  directly  treating  the 
subject  of  the  bow  shock  and  magnetopause  motion  induced  through  collision  with  an 
interplanetary  shock  wave.  The  first  article,  by  Shen  and  Diyer  (1972),  discusses  the 
interaction  of  an  interplanetary  shock  wave  with  the  earth’s  bow  shock  and  provides  a 
technique  to  solve  for  the  propagation  speeds  of  the  resulting  shock  waves.  The  second 
article,  by  Shen  (1973)  alone,  is  a  qualitative  study  of  the  interaction  of  a  shock  wave 
with  the  magnetopause.  The  third  article  is  a  detailed  work  by  Grib,  et.  al.  (1979)  which 
develops  analytic  solutions  for  shock-bow  shock  and  shock-magnetopause  interactions  in 
both  the  gas-dynamic  and  MHD  cases.  Grib’s  solutions  provide  a  means  to  analytically 
determine  position  and  velocity  for  the  displaced  bow  shock  and  magnetopause  as  a 
fimction  of  the  Mach  number  of  the  incident  interplanetary  shock  wave;  they  can  also 
serve  as  a  validation  tool  to  asses  the  ability  of  the  numerical  model  developed  in  this 
thesis  to  predict  bow  shock  and  magnetopause  motion. 

3.1  Shen  and  Dryer 

Shen  and  Dryer  (1972)  formulated  the  interaction  as  a  one  dimensional  problem 
in  which  the  interplanetary  shock,  a  fast  magnetosonic  wave  propagating  through  the 
solar  wind  from  the  sim  to  the  earth,  collides  with  a  standing  bow  shock  (1972:4628-29). 
Echoing  earlier  treatments  by  Courant  and  Friedrichs,  and  Landau  and  Lifshitz,  Shen  and 
Dryer  discuss  the  formation  of  two  new  shock  waves,  S3  and  S4,  with  a  contact 
discontinuity,  C],  between  them  generated  by  the  ISW-bow  shock  interaction.  Applying 
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continuity  of  flow  velocity  and  pressure  (Equations  (24)  and  (26),  with  B  =  0)  across  Ci 
for  the  gas-dynamic  case,  they  obtain  the  following  equations; 


and. 


2yM/-(Y-1)][2yMi^-(Y-1)] 

2yM3'*-(y-1)][2yM22-(y-1)]  ■  ^ 


l-M/ 
M,(y  +  1) 


1-M,^ 
M,(y  +  1) 


where. 


G(M,y) 


fpYM^-(Y-l)l  r2+M"(Y-l)11*'' 

U  Y  +  1  Jl  (Y+1)M^  Jj 


is  the  square  root  of  the  ratio  of  the  shock  strength  and  compressibility,  "'/(y/ri),  for 
specified  Mach  number,  M,  and  isentropic  exponent,  y  (Shen,  1972:4631).  Simultaneous 
solutions  of  Equations  (27)  and  (28)  yield  the  Mach  numbers,  Mj  and  M4,  of  the  two 
generated  shocks  in  terms  of  the  known  bow  shock  and  ISW  Mach  numbers,  Mj  and  M2. 

To  more  realistically  describe  the  ISW-bow  shock  interaction,  Shen  and  Dryer 
include  magnetic  effects  by  applying  the  MHD  form  of  Equations  (24)  and  (26)  (B  ^  0); 
continuity  of  flow  velocity  and  total  pressure  across  C]  are  then  expressed  as. 


A,p(M„A,)  +  A,P(M„A,)[^(M.,A,)f 

-  A2P(M2,A2)  +  A3P(M3,A3)[^(M2,A2)f  =  0 
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and. 


2  Mo"  2  A," 

5(M„a,)5(M„AJ  ^"(M„A,)^"(M4,AJ 

2M,"  2A," 


where. 


A  =  AlfVen  number  =  the  ratio  of  the  wave  speed,  V,  to  the 
Alfven  speed,  Va, 

a"+m"-a"  m" 


P(M,A)  = 


3A"M" 


3A"M" 


^(M,A)  = 


2M"  +2A"  +M"A"’ 


with  A3  and  At  given  by, 

^  _  M.  f5(M3,A3)]'^' 

^  A,  4(M3,A3) 

A  _  M,  [5(M„A,)f 
^  A,  ^(M„A,) 


(31) 

(32) 


Shen  and  Dryer  define  ^(M,A)  as  the  MHD  compressibility  and  6(M,A)  as  the  MHD 
shock  strength  (Shen,  1972:4633).  The  simultaneous  solutions  of  Equations  (29)  -  (32) 
yield  acoustic  Mach  numbers  M3  and  Mi,  and  the  Afv6n  numbers,  A3  and  At,  for  shock 
waves  S3  and  S4. 
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With  the  Mach  and  AlfVen  numbers  of  the  generated  shocks  determined,  then 
relationships  given  by  Shen  and  Dryer  in  Table  2  can  be  used  to  determine  the 
propagation  speeds  of  S3  and  S4  through  the  magnetosheath. 


Shock 


S, 

52 

53 

54 


Table  2.  Mach  number  relations  for  shock-bow  shock  interaction  _  _ 

Acoustic  AlfVen  Magnetosonic 

Mach  Number  _ Numbw _ MachN^ber 

'  “m,  =  (uo-V,)7ao  Ai  =  (uo-V,)/b7  '  M*r=  (uo-V,)/co  “ 

Mi  =(uo-V2)/ao  A2  =(uo-V2)/bo  M*2  =(uo-Vi)/co 

M3  =  (U2-V3)/a2  A3  =  (U2-V3)/b2  M*3  =  (U2-Vr)/C2 

M4  =  (ui-V4)/ai  A4  =  (ui-V4)/bi  M*4  =  (ui-V4)/c| 

from;  Shen  and  Dryer  (1972:4633) 


where  the  form  is  M*  =  (uj  -  Vjt)/ai,  and, 

u,  =  gas  flow  velocity  in  Region  /  of  Figure  4,  /  =  0, 1, 2, 

Yk  =  shock  propagation  speed  of  the  Ath  shock,  A:  =  1, 2, 3, 4, 

a,  =  speed  of  soimd  Eq  (1)  in  Region  /, 

b,  =  AlfVen  speed  in  Region  i, 

c,  =  fast  magnetosonic  sound  speed  Cf  Eq  (7)  in  Region  /. 

These  relationships  hold  in  the  absolute,  or  earth-centered,  reference  frame.  Note  that 
the  absolute  frame  is  also  the  bow  shock  reference  frame  because  Vi,  the  propagation 
speed  of  the  bow  shock  before  collision  with  the  ISW,  is  taken  to  be  equal  to  zero. 


3.2  Shen 

In  his  second  paper,  Shen  (1973)  considers  the  interaction  of  S4  with  a  stationary 
magnetopause  characterized  as  a  tangential  discontinuity  (1973:55).  He  again  proceeds 
along  the  same  lines  as  Courant  and  Friedrichs,  and  Landau  and  Lifshitz,  by  qualitatively 
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describing  the  formation  of  a  new  shock  wave,  Se,  and  a  reflected  rarefaction  wave,  R5, 
resulting  from  the  S4-magnetopause  collision  (see  Figure  5).  Invoking  the  continuity  of 
plasma  flow  velocity  and  pressure  across  the  new  contact  discontinuity,  0^,  also  formed 
in  this  collision,  Shen  states  that  a  system  of  four  equations  in  four  unknowns  can  be 
solved  to  yield  the  Alfv^n  number  for  R5  and  the  Mach  number  for  Se;  he  does  not, 
however  provide  the  equations  to  complete  the  solution.  He  does  provide  a  useful  table 
of  relationships  for  the  AlfVenic  and  acoustic  Mach  numbers  presented  in  Table  3. 

Table  3.  Mach  number  relations  for  shock-magnetopause  interaction 


Shock 

Ordinary 

Mach  Number 

Alfven 

Number 

Magnetosonic 
Mach  Number 

S4 

M4  =  (V4-u,)/a, 

A4  =  (V4-ui)/bi 

M*i  =  (V4-u,)/c, 

R5 

M5  =  (V5-U4)/a4 

As  =  (V5-U4)/b4 

M*5  =  (V5-U4)/C4 

S6 

M6  =(V6-u„.)/a„ 

^6  ('^6"trin)^m 

M*6 

from  Shen  (  1973:54) 

where  the  subscripts  now  refer  to  the  regions  in  Figure  6,  and, 

u,  =  gas  flow  velocity  in  Region  i  of  Figure  6,  /  =  1, 4,  and  m, 

V*  =  wave  propagation  speed  of  Ath  wave.  A:  =  4,  5,  and  6. 

Two  additional  points  discussed  by  Shen  deserve  mention  here.  The  first  is  that 
observations  indicate  density  “is  sharply  cut  off  at  the  boundary  of  the  magnetospheric 
cavity”  (1973:56),  further  supporting  that  the  case  of  interest  when  S4  interacts  with  the 
magnetopause  is  the  generation  of  a  reflected  rarefaction  wave  (recall  Courant  and 
Friedrichs).  The  second  point  is  that  because  of  lower  density  and  a  small  plasma  (3 
inside  the  magnetosphere,  the  propagation  speed  of  Sg  inside  the  magnetosphere  is 
substantially  higher  than  the  velocity  of  the  shock  wave  S4  (1973:59). 
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3.3  Grib 


The  third  and  final  paper  reviewed  is  an  excellent  treatise  by  S.  A.  Grib,  et.  al., 
(1979)  on  the  entire  problem  from  the  ISW  collision  with  the  bow  shock  to 
magnetopause  retrograde  motion  toward  its  initial  position.  This  paper  is  mathematically 
rich  and  most  of  the  equations  derived  therein  are  important  to  the  development  of  this 
thesis;  therefore  many  of  their  equations  will  be  included  in  this  review  for  reference  in 
the  next  chapter.  To  ensure  that  proper  credit  is  given  to  other’s  work  and  yet  avoid 
repeated  reference  to  the  same  source,  subsequent  citations  within  this  chapter  from  Grib, 
et.  al.  (1979)  will  reference  only  the  page  number  unless  otherwise  stated. 

As  Shen  and  Dryer  before  them,  Grib,  et.  al.  assume  a  one-dimensional  problem. 
They  also  state  that  only  the  component  of  the  magnetic  field  perpendicular  to  the  flow 
velocity — i.e.  in  the  z  direction — is  significant  to  the  problem  (5908),  thus  the  fast 
magnetosonic  and  acoustic  waves  are  the  only  waves  of  concern.  As  preparatory 
material  they  also  give  two  useful  equations,  based  on  the  Riemann  invariants,  for  plasma 
flow  speed  behind  a  rarefaction  wave  in  the  gas-dynamic  (Eq  (33))  and  MHD  (Eq  (34)) 
cases: 


-2  _  2 
u+ - -a  =  + - -ao  =  const 

Y-1  Y-1 

2  2 
u±— aJ,(a)  =  ±— ao  J,(a) 

where  the  subscript  0  refers  to  upstream,  “unshocked”  values,  and 
u  s  gas  flow  speed. 


(33) 


(34) 
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2  2  2  2 

oc=VA/a  =B/  2iiopa  =  the  ratio  of  AlfVen  speed  to  sound  speed, 

1 

Jg(a)  s  JVl+aj^^^dx . 

0 

Proceeding  along  the  same  lines  as  Shen  and  Dryer,  Grib  describes  the  creation  of 
the  two  new  shocks,  S3  and  84^  and  the  contact  discontinuity  between  them,  resulting 
from  the  ISW-bow  shock  collision  as  illustrated  in  Figure  4.  Given  a  Mach  number  for 
the  ISW  and  ambient  solar  wind  conditions,  then  the  flow  parameters  for  plasmas  in 
Region  1  and  Region  2  are  easily  determined  from  the  shock  jump  equations.  Assuming 
that  the  bow  shock  is  stationary  in  the  earth  reference  frame  determines  that  the  upstream 
flow  velocity  into  the  bow  shock  is  just  the  ambient  solar  wind  speed;  Region  1  values 
are  constant  for  a  given  ambient  solar  wind  and  do  not  change  until  disturbed  by  passage 
of  S4.  Parameters  behind  the  incoming  ISW  depend  on  the  Mach  number  M2  and  on  the 
conditions  of  the  upstream  solar  wind,  but  for  any  given  solar  wind  conditions  parameters 
in  Region  2  depend  only  on  M2. 

Applying  equality  of  pressure  and  velocity  across  the  contact  discontinuity,  and 
assuming  that  U2,  P2,  Uj,  and  pi  are  known  variables,  Grib  arrives  at  a  system  of  two 
equations  in  M3  and  M4  (for  the  gas-dynamic  case): 


2a2  Mj^-l  2a,  M/-1 
y-1  M3  M, 


(35) 

(36) 
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Here,  the  left  hand  sides  of  Equations  (35)  and  (36)  are  the  flow  speed  and  pressure  in 
Region  3;  the  right  hand  sides  are  the  flow  speed  and  pressure  in  Region  4.  Solving  these 
equations  simultaneously  result  in  solutions  of  Mj  and  H,  in  terms  of  M2.  Grib  also 
defines  the  Mach  number  of  the  plasma  flow  in  Region  2  as  M*  =  U2  /  a2.  The  numerical 
solutions  of  Equations  (35)  and  (36).are  plotted  in  Figure  6,  for  Mach  Numbers  M3,  Mt, 
and  M*  in  terms  of  M2  for  the  gas-dynamic  case.  ‘ 


M2 

Figure  6.  Mach  numbers  for  gas-dynamic  ISW-bow  shock  interaction. 

M3  is  the  Mach  number  for  the  reflected  shock,  M4  is  the  Mach  number 
of  the  transmitted  shock,  and  M*  is  the  Mach  number  of  the  fluid  flow 
downstream  of  the  ISW. 

Knowing  M4  and  aj,  Grib  develops  an  equation  to  determine  the  time,  t2  -  ti, 
required  for  S4  to  cross  the  magnetosheath  and  reach  the  magnetopause.  In  the  earth- 
centered  reference  frame,  the  new  bow  shock,  S3,  travels  toward  the  earth  with  a  velocity 


34 


a2(M*  -  M3)  (5909);  until  S4  reaches  the  magnetopause  at  time  ta,  the  thickness  of  the 
magnetosheath  decreases  as  S3  moves  earthward.  Defining  the  initial  thickness  of  the 
magnetosheath  as  6,  then  the  change  in  thickness  at  some  At  =  t2  -  ti  is 

At 

|5, -5|  =  Ja2(M*-M3)dt  (37) 

0 

where, 

81  is  the  new  magnetosheath  thickness, 
and  At,  the  time  it  takes  S4  to  traverse  the  magnetosheath,  given  by, 

_ dx _  6  f  u,  5 

|(u,(l-8/x)  +  a,M4)  u,  ^a,M4J  ~  a,M4  ■ 

(5910) 

Eq  (38)  is  based  on  the  assumption  that  the  flow  velocity  in  the  imdisturbed 
magnetosheath  reduces  linearly  to  zero  at  the  magnetopause.  Eq  (37)  allows  the  position 
of  the  moving  bow  shock  to  be  determined  as  the  change  in  magnetosheath  thickness  in  a 
time  At. 

Having  determined  the  Mach  numbers  M3  and  M4,  and,  consequently,  all  of  the 
plasma  parameters  in  Regions  3  and  4  by  application  of  the  jump  conditions,  Grib  turns 
to  a  solution  of  the  S4-magnetopause  interaction.  Before  S4  arrives  at  the  magnetopause, 
there  must  be  a  balance  of  pressures  between  the  magnetosheath  and  the 
magnetosphere  a  boundary  condition  for  a  tangential  discontinuity  (recall  Parks, 

1991 .329).  However,  observations  show  that  density  inside  the  magnetosphere  is 
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significantly  lower  than  in  the  magnetosheath  so  that,  from  the  ideal  gas  law,  pressure  in 
the  magnetosphere,  is  much  lower  than  pi  and  the  gas  dynamic  approach  fails. 
Because  dynamic  pressure  alone  in  the  magnetopause  is  insufficient  to  balance  the 
pressure  from  the  magnetosheath,  Grib  adopts  a  quasi-gas-dynamic  approach  in  which 
the  dynamic  pressure  of  the  solar  wind  is  balanced  by  the  sum  of  dynamic  and  magnetic 
pressures  within  the  magnetosphere.  By  including  the  magnetic  pressure  of  the  earth’s 
relatively  strong  planetary  field,  pressure  balance  is  achieved  across  the  magnetopause 
even  though  p^  is  less  than  pi  (591 1). 

Collision  of  S4  with  the  magnetopause  produces  a  new  shock,  Se,  transmitted  into 
the  magnetopause,  a  reflected  rarefaction  wave,  R5,  propagating  sunward  toward  the 
incoming  bow  shock,  and  a  new  contact  discontinuity,  which  in  effect  becomes  the 
new  magnetopause  moving  with  the  flow  velocity  in  Region  5  and  6  toward  the  earth  (see 
Figure  5).  The  rarefaction  wave  is  not  a  shock  and  therefore  propagates  as  an  acoustic 
wave  (5908)  into  Region  4.  In  the  absolute  reference  frame  the  rarefaction  wave  travels 
at  velocity  a4  -  U4  so  that  R5  will  travel  toward  the  advancing  bow  shock  only  if  a4  >  U4. 
For  solar  wind  parameters  listed  in  Table  1,  a4  >  U4  if  M2  <  5  (5911).  It  will  be  shown  at 
the  end  of  this  chapter  that  the  interaction  of  rarefaction  waves  is  a  mechanism  that 
reverses  the  earthward  motion  of  the  bow  shock  and  magnetopause.  Thus,  when  M2  >  5, 
the  rarefaction  wave  can  not  propagate  sunward  to  interact  with  the  advancing  bow  shock 
and  eventually  reverse  its. 

At  the  new  magnetopause.  Cm,  continuity  of  total  pressure  and  flow  velocity,  and 
the  Riemann  invariants  associated  with  the  rarefaction  wave  (Eq  (33))  results  in  a  set  of 
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two  equations  in  terms  of  ps,  the  density  in  Region  5,  and  M^,  the  Mach  number  of  the 
shock  wave  propagated  into  the  magnetosphere; 


and. 


U4  + 


2 

Y-1 


(39) 


P. 


2  Y  M4^  + 1  -  y 

P5 

■(Y-1)M42+2' 

y  f 

Y  +  1 

-Pi. 

(Y  +  1)M4 

V 

1-C 

1+1/p 


(40) 


where  C  =  C  (M^,  p„)  (Eq  (19)),  and. 


£1 

^pj 


(r-i)/2 


(41) 


P,  =  P,(^J.  (42) 

(5911) 

Knowing  the  Mach  number  of  an  approaching  interplanetary  shock  wave,  M2,  and 
the  conditions  of  the  solar  wind.  Equations  (35),  (36),  and  (39)  -  (42)  permit  analytic  gas- 
dynamic  solutions  for  Mj,  M4,  and  M^  in  terms  of  M2.  These  Mach  numbers  are  then 
used  to  determine  bow  shock  and  magnetopause  motions,  as  well  as  plasma  parameters 
in  each  region  of  Figures  4  and  5. 

Just  as  Shen  and  Dryer  before  them,  Grib,  et.  al.  repeat  their  analysis  for  the  MHD 
case  by  including  the  effects  of  magnetic  fields  on  shock  wave  propagation  in  a  plasma. 
Here  the  full  MHD  forms  of  the  Rankine-Hugoniot  equations  and  shock  jump  equations 
must  be  used.  The  physical  arguments  leading  to  the  formation  of  two  new  shocks  upon 
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collision  of  the  ISW  with  the  bow  shock,  and  the  rarefaction  wave-transmitted  shock  pair 
resulting  from  the  S4-magnetopause  interaction,  do  not  change  in  the  MHD  approach  so 
that  Figures  4  and  5  still  represent  these  interactions.  Different  forms  of  the 
compressibility,  and  the  shock  strength,  y,  lead  to  different  jump  equations  and  to 
different  relationships  for  pressure  and  velocity  continuity  across  a  contact  discontinuity. 
Grib  expresses  the  relationships  for  continuity  of  pressure  and  velocity  at  the  contact 
discontinuity  Ci  as 


U2-a2M3(l-Cj)  =  u,+a,M4(l-C4) 


(43) 


P2 


f  l_r 
1+— !  yM," 


V  1+1/P2 


I-C4 
1  +  1/p, 


(44) 


where, 

p*  =  the  total  pressure  (Eq  ( 1 )), 

=  compressibility  ratio  in  terms  of  the  Mach  number  M  and  p  upstream  of  the 
shock  (Eq  (19)), 

P  s  the  ratio  of  magnetic  to  dynamic  pressure  =  B^/(2  pop)  (5913). 

Similar  to  the  gas-dynamic  treatment.  Equations  (43)  and  (44)  are  solved  simultaneously 
to  obtain  M3  and  M4  for  the  MHD  case  and  the  results  are  presented  in  Figure  7. 
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Figure  7.  Mach  numbers  for  MHD  ISW-bow  shock  interaction. 


The  solution  for  M4  and  the  plasma  conditions  behind  the  shock  S4,  allow  Grib  to 
make  the  next  step  in  the  analysis;  an  examination  of  the  flow  velocities  and  total 
pressures  in  Regions  5  and  6  that  will  enable  a  solution  for  and  the  plasma  parameters 
behind  the  transmitted  shock.  At  the  new  magnetopause.  Cm,  equality  of  flow  speed  is 
given  by 


(l-C)a,M,+;^[a,  J.(o,)-a,J,(a,)]  =  (l-Qa,M.  (45 ) 


where. 


ttj  = 


^0P4C4  VpJ 
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V  p  ; 


and,  continuity  of  pressure  gives, 


B. 


T  f .  N2 


2^io 


£i 

VPi. 


=  Pm 


1-C6 


(46) 

(:5913) 


Solving  Equations  (45)  and  (46)  simultaneously  results  in  a  determination  of  the  two 
unknowns,  ps  and  Mg.  Equations  (43)  -  (46)  form  a  set  of  simple  analytical  relationships 
from  which  the  MHD  Mach  numbers  for  all  three  generated  shocks  are  found  in  terms  of 
M2.  With  these  Mach  numbers,  along  with  the  ambient  conditions  in  the  solar  wind  and 
magnetopause,  plasma  parameters  in  each  region  of  Figures  4  and  5  are  determined  for 
the  MHD  case.  The  new  magnetopause  will  move  with  the  flow  at  speed  Ug  after 
collision  with  S4;  the  new  bow  shock,  S3,  moves  with  a  speed  U2  -  (1  +  M3)a2  (Grib,5913). 

Having  determined  the  conditions  of  the  magnetosheath  up  to  a  time  just  after  the 
creation  of  Sg  and  the  rarefaction  wave  R5,  Grib  et.  al.  take  their  examination  one  step 
farther  by  determining  when  and  how  the  magnetopause  will  reverse  its  earthward 
movement.  They  describe  a  series  of  successive  rarefaction  waves  reflecting  off  of  the 
advancing  bow  shock  rear  and  the  magnetopause  as  depicted  in  Figure  8.  The  first 
rarefaction  wave,  R5,  decreases  the  particle  density  and  pressure  in  the  magnetosheath  as 
it  travels  against  the  flow  toward  the  bow  shock.  It  will  reach  the  bow  shock  only  if  its 
speed  is  greater  than  the  flow  velocity  in  the  magnetosheath,  requiring  M2  >  6  for  Bq  = 

3.5  X  10'^  Teslas  (nT)  (:5915).  At  the  rear  of  the  bow  shock,  a  new  rarefaction  wave,  or 
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rarefaction  fan,  labeled  Rg  reflects  off  of  the  bow  shock  and  moves  toward  the 
magnetopause;  the  bow  shock  speed  increases  slightly  as  it  advances  into  the  region  of 
lower  density  caused  by  R5.  Traversing  the  magnetosheath,  Rg  encounters  the 
magnetopause  where  it  decreases  the  flow  pressure  just  behind  Cm,  allowing  the 
magnetopause  to  slow  and  to  reverse  direction  moving  back  against  the  flow.  The 
reversed  magnetopause  now  “pushes”  on  the  gas  in  front  of  it  much  like  a  piston  in  a 
shock  tube,  creating  multiple  small  perturbations  in  pressure  that  steepen  into  a  new 
shock  which  interacts  with  the  bow  shock  to  produce  a  new  bow  shock  Sjo  (:5914-15). 


Cm 


Figures.  Rarefaction  wave  interactions  (Grib,  1979:5911) 

In  their  paper,  Grib  and  his  colleagues  arrive  at  an  equation  of  magnetopause 
motion  in  the  form. 
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X  =  (U6-C6)t  +  (C6+C5)ty 


-C5T 


(47) 


where, 

X  s  time  between  arrival  of  the  shock  S4  and  the  arrival  of  Rs  =  h-  tj, 

c  =  fast  magnetosonic  speed  in  Regions  5  and  6,  respectively  (:5915). 

Consistent  with  Grib’s  explanation  of  magnetopause  reversal,  Figure  8  is  drawn  to  depict 
the  magnetopause  turning  to  coincide  with  the  intersection  of  Rg  which  occurs  at  time  x 
after  the  initial  shock  impact.  Grib,  et.  al.  provide  estimated  limits  for  x  based  upon  the 
magnetic  field  of  the  solar  wind,  Bq.  For  Bq  =  0  x  =  1-3  minutes;  for  Bq  =  3.5  nT  x  =  3-5 
minutes;  and  for  Bq  =  7.0  nT  x  =  7  minutes  (Grib,  1979:5914). 

3.4  Literature  Review  Summary 

In  this  literature  review  three  articles  were  discussed:  Shen  and  Dryer’s  treatment 
of  the  ISW-bow  shock  interaction,  Shen’s  paper  on  the  S4-magnetopause  interaction,  and 
the  detailed  analysis  of  the  complete  ISW-magnetosheath  interaction  presented  by  Grib, 
et.  al.  From  Shen  and  Dryer  comes  a  technique  for  solving  sets  of  gas-djmamic  and 
MHD  equations  for  the  unknown  Mach  numbers  of  the  shock  waves  generated  by  the 
ISW-magnetosheath  interaction;  i.e.  application  of  the  continuity  of  plasma  flow  velocity 
and  pressures  across  a  contact  discontinuity  and  simultaneous  solution  of  the  resulting 
equations  for  the  Mach  numbers.  Shen’s  second  paper  is  a  qualitative  discussion  of  the 
shock-discontinuity  interaction.  In  the  article  by  Grib,  et.  al.,  the  technique  of  Shen  and 
Dryer  is  applied  to  both  the  ISW-bow  shock  and  the  S4-magnetopause  interactions  to 
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yield  sets  of  gas-dynamic  and  MHD  equations  which  are  then  solved  to  determine  the 
unknown  Mach  numbers  of  the  generated  shock  waves.  Applying  the  Mach  number 
solutions  to  the  shock  jump  equations  allows  plasma  parameters  in  each  region  of  the 
modified  magnetosheath  to  also  be  determined,  leading  to  solutions  for  bow  shock  and 
magnetopause  motion.  Extending  their  treatment  beyond  the  initial  motions  of  the  bow 
shock  and  magnetopause,  Grib,  et.  al.  arrive  at  an  equation  of  magnetopause  motion 
which  describes  the  reversal  of  the  magnetopause  movement  resulting  from  the  repeated 
interaction  of  rarefaction  waves  between  the  bow  shock  and  magnetopause. 
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IV.  Methodology 


The  earlier  works  of  Shen ,  Dryer,  and  Grib  lay  a  qualitative  and  analytic 
foundation  toward  understanding  bow  shock  and  magnetopause  displacement  resulting 
from  collision  with  an  interplanetary  shock  wave.  To  advance  the  understanding  of  the 
physical  processes  associated  with  bow  shock  and  magnetopause  motion,  and  to 
approximate  to  a  reasonable  degree  of  confidence  the  displacements  and  velocities 
resulting  from  these  processes,  a  one-dimensional  numerical  model  is  developed  to 
simulate  the  ISW-magnetosheath  interactions.  This  chapter  is  a  presentation  of  the 
development,  validation,  and  execution  of  this  numerical  model. 

Development  of  the  numerical  model  begins  with  the  fluid  conservation 
equations,  from  which  a  set  of  time  and  spatially  dependent  differential  equations  for  the 
conservation  variables — i.e.  mass  density,  momentum,  magnetic  field,  and  total  energy — 
are  derived.  Next,  the  numerical  algorithm  used  for  a  one-dimensional,  explicit,  time 
dependent  solution  of  the  differential  equations  is  discussed.  Procedures  used  to  extract 
pressure,  density,  temperature,  magnetic  field,  and  flow  velocity  from  conservation 
variables  updated  by  the  model,  as  well  as  the  method  for  determining  shock  front  and 
contact  discontinuity  locations  at  each  time  step  are  also  presented. 

After  detailing  the  development  of  the  numerical  model,  the  next  discussion 
examines  the  procedures  used  to  validate  the  numerical  model.  To  verify  that  the 
algorithm  generates  reasonable  and  expected  results,  several  test  cases  are  executed  and 
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the  results  are  analyzed  and  compared  to  shock  propagation  speeds  derived  from  known 
analytic  solutions. 

The  final  section  in  Chapter  IV  describes  model  input,  initialization  procedures, 
and  data  processing. 

4.1  Model  Development 

There  are  several  advantages  of  adopting  a  numerical  approach  to  determine  bow 
shock  and  magnetopause  displacement  when  analytic  solutions  already  exist.  One 
advantage  is  that  a  numerical  simulation  allows  observation  of  the  time  evolution  of  the 
shock-shock  and  shock-discontinuity  interactions.  Another  advantage  is  that  a  numerical 
simulation  quickly  provides  a  reasonable  approximation  of  not  only  bow  shock  and 
magnetopause  displacement,  but  it  also  yields  information  about  plasma  conditions — i.e. 
pressure,  density,  temperature,  magnetic  field,  and  flow  velocity— for  each  region  of  the 
interacting  system  at  each  time  step.  As  a  consequence  of  this  ability  to  determine 
plasma  parameters  in  different  regions,  a  numerical  simulation  can  simultaneously  track 
the  progression  of  several  different  shock  waves,  contact  discontinuities,  and  rarefaction 
waves.  A  third  advantage  is  that  the  complicated  interactions  between  rarefaction  waves 
and  the  bow  shock  and  magnetopause  described  by  Grib,  et.  al.  can  theoretically  be 
simulated  by  numerical  code.  To  arrive  at  their  equation  of  magnetopause  motion  (Eq 
(47)),  Grib  and  his  colleagues  made  evaluations  of  rarefaction  wave  propagation  based 
on  average  velocities  of  the  waves  through  different  plasma  regions  (Grib,  1979:5912, 
5914);  by  tracking  changes  in  plasma  parameters  characteristic  of  a  rarefaction  wave, 
numerical  simulations  would  provide  better  estimates  of  rarefaction  wave  position  and 
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velocity. 


4.1.1.  Equation  Development.  Development  of  a  numerical  model  began  with 
the  vector  MHD  fluid  equations  given  by  Tdth  (1996:83)  reduced  to  the  one  dimensional 
case: 


ap  d  /  \ 

^  +  -^(pif)=0 

dt  ax'  ' 

(48) 

>4(uB).0 

(49) 

±(puh£(pua)  =  -^(p) 
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where, 
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p*  =  total  pressure  (Eq  ( 1 )). 

Substituting  e  and  p*  into  the  fluid  equations.  Equations  (48)  -  (51)  become. 
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Equations  (52)  -  (55)  form  a  set  of  time  and  spatially  dependent  partial  differential 
equations  in  the  conservation  variables.  This  set  of  equations  can  be  written  in  a 
convenient  matrix  form, 


where. 
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It  is  often  desirable  to  non-dimensionalize  a  set  of  equations  used  in  a  numerical 
model  to  simplify  calculations.  Equations  (57)  -  (58)  are  non-dimensionalized  by 
applying  the  parameter  transformations  listed  below: 
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X  = 
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where, 

L  =  scale  length  of  the  model, 

Vfc  =  thermal  velocity  of  the  solar  wind  =  (k  To/m)'^,  and  the  superscript  0  refers 
to  the  ambient  solar  wind  values  of  the  indicated  parameters.  The  non-dimensionalized 
Q  and  F  matrices  are. 
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where  the  primed  notation  signifies  dimensionless  parameters  and  Po  is  defined  as  the 
ratio  of  magnetic  to  dynamic  pressure,  Bo^  /  (2  po  Po),  in  the  solar  wind.  Note  that  the 
factor  2  multiplying  pressure  and  magnetic  field  terms  in  Equations  (59)  and  (60)  is  due 
to  an  ideal  gas  law  in  the  form  p  =  2  n  k  T  which  accoimts  for  a  quasi-neutral  plasma  in 
which  the  number  densities  and  temperatures  of  ions  and  electrons  are  assumed  equal. 
The  ideal  gas  law  becomes  p  «  2  p  k  T/m,  where  p  is  the  ion  mass  density  and  m  is  the 
ion  mass,  when  the  mass  of  the  ions  are  much  greater  than  the  mass  of  the  electrons. 

The  function  of  the  algorithm  is  to  numerically  update  Q  in  time.  To  determine 
shock  and  discontinuity  locations,  however,  plasma  pressure,  density,  temperature, 
magnetic  field,  and  flow  velocity  must  be  known  at  each  bin  in  the  model  at  successive 
time  steps.  These  quantities  are  calculated  from  values  of  Q  by  the  following  equations; 


P'  = 


P'  =  Q(i) 

(61) 

B'  =  Q(2) 

(62) 
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Q(i) 
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The  numbers  in  parentheses  refer  to  elements  of  the  Q  matrix. 
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Equations  (48)  -  (65)  are  valid  for  both  the  gas-dynamic  and  MHD  fluid  approximations. 
If,  in  these  equations  B  =  0,  then  magnetic  terms  drop  out  and  the  gas-dynamic  equations 
is  recovered. 

4.1.2.  Algorithm  Development  With  the  derivation  of  the  governing  fluid 
equations  complete,  the  next  step  was  the  development  of  a  suitable  numeric  algorithm  to 
advance  solutions  of  the  equations  in  time.  The  core  algorithm  selected  was  a 
MacCormack  predictor-corrector  scheme  with  Flux  Corrected  Transport  (FCT)  taken 
from  Fletcher’s  book.  Computational  Techniques  for  Fluid  Dynamics:  Vol  II,  Techniques 
for  Different  Flow  Categories  (1991:154-156).  This  finite  difference  technique  is 
effective  in  capturing  shock  locations  within  a  steady  inviscid  supersonic  flow 
(1991:147);  it  is  also  fairly  computationally  simple  and  can  be  implemented  on  a  desktop 
computer.  FCT  is  a  technique  incorporated  into  the  algorithm  to  reduce  numerical 
oscillations  at  locations  where  steep  gradients,  such  as  shock  waves  and  discontinuities, 
cause  numerical  instabilities.  The  technique  adds  higher  order  diffusive  terms  to  broaden 
a  shock  or  discontinuity  over  a  greater  number  of  spatial  bins,  thus  decreasing  the 
gradient  and  the  associated  instabilities.  Higher  order  anti-diffusive  terms  are  then 
subtracted  out  of  the  diffused  solution  to  recover  a  physically  meaningful  result.  The  art 
in  the  technique  is  to  add  diffusive  terms  and  subtract  anti-diffusive  terms  in  just  the  right 
amounts  (Toth,  1996:82). 

The  MacCormack  algorithm  solves  the  conservation  equations  of  mass,  magnetic 
flux,  momentum,  and  energy  in  the  form. 
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(56) 


dQ  ^ 
dt  dx 

across  each  bin  advancing  Q  in  time  by  two  stages  (Fletcher,  1991: 148).  Here  Q  and  F 
are  the  matrices  derived  in  Equations  (59)  and  (60).  The  first,  or  predictor  stage,  obtains 
an  intermediate  solution 


which  is  then  put  into  the  second,  or  corrector  stage 

q,-  =  o.5(q-+q;)-^(f;-f;,)  (6?) 

where, 

superscript  n  designates  the  value  at  time  tn  =  ^^nAtj , 

subscript)  is  the  spatial  bin  index. 

Ax  =  spatial  bin  width  in  dimensionless  units. 

At  =  optimized,  dimensionless  time  step  for  time  t„, 

Q*  is  the  intermediate  solution  of  the  predictor  stage, 

F*  is  the  intermediate  value  of  F  as  a  fimction  of  Q*.  (Fletcher,  1991:149). 

The  time  step.  At,  must  be  constrained  to  ensure  that  gas  flow  does  not  advance 
more  than  one  bin  in  a  single  step,  otherwise  the  model  may  become  unstable.  The 
desired  time  step  is  one  in  which  [  |  uj  |  ]max  ^  Ax,  where  [  |  Uj  |  ]max  is  the  absolute 
value  of  the  speed  taken  from  the  bin  with  the  largest  flow  velocity.  To  guarantee  that 
flow  advances  at  most  only  one  bin  in  a  time  step,  the  right  hand  side  of  the  inequality  is 
multiplied  by  a  factor  less  than  one  called  the  Courant-Friedrich-Levy  (CFL)  factor  (Toth, 
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1996:84),  which  in  the  present  code  is  a  constant  equal  to  0.25.  From  Toth,  (1996:84), 
the  time  required  to  step  the  model  so  that  changes  in  mass,  momentum,  and  energy 
fluxes  propagate  no  more  than  one  bin  is, 

it  s  CFLjj— •  («*) 

I  "i 

The  denominator  in  Eq  (68)  is  the  wave  propagation  speed  in  the  absolute  reference 
frame  taken  from  the  bin  having  the  largest  absolute  value  of  |  uj  |  +  Cq,i.  Here  Cq,  the 
speed  of  the  fast  mode  MHD  wave,  is 
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(Toth,  1996:84). 


In  dimensionless  form  Eq  (69)  becomes. 
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Equations  (69)  -  (70)  are  valid  for  both  the  gas-dynamic  and  MHD  cases.  If  B  =  0,  then 
magnetic  terms  drop  out  and  Cq  becomes  just  the  acoustic  wave  speed. 

Following  stages  one  and  two,  Q  is  sent  to  the  FCT  subroutine  where  a  diffusive 
flux,  DF,  is  added  to  reduce  the  gradients  in  parameter  jumps,  (Eq  (68)) 


Qj 


=  Q^+DF3-DFj., 


(71) 
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and  anti-diffusive  fluxes,  ADFs,  are  subtracted  (Eq  (69))  to  recover  the  desired  solution. 

Qj=Q3-hADF3.,-ADFj  (72) 

The  details  of  how  DF  and  ADF  are  calculated  provided  in  Fletcher’s  book. 

The  last  step  in  the  MacCormack  algorithm  is  to  extract  from  the  new  gas 
parameters  in  each  bin,  which  is  accomplished  using  Equations  (61)  -  (65)  derived  in 
Section  4.1.1. 

An  enhancement  to  the  basic  code  was  the  inclusion  of  a  new  procedure  to  locate 
multiple  structures — shocks  and  discontinuities — from  the  data  extracted  from  A 

density  jump,  or  large  density  gradient,  is  a  feature  common  to  both  shock  waves  and 
discontinuities  and  is  relatively  easy  to  extract  from  algorithm  data.  A  new  subroutine 
was  written  to  locate  bins  where  the  density  gradient,  Ap/Ax,  exceeded  a  threshold  value 
of  100,  corresponding  to  a  density  change  of  approximately  four  percent  between 
adjacent  bins.  The  routine  begins  at  the  largest  bin  number  and  progresses  sequentially 
backward  (from  right  to  left)  testing  density  gradients  at  each  bin.  When  the  gradient 
exceeds  the  threshold  value,  the  routine  enters  a  second  loop,  beginning  at  the  current 
bin,  and  continues  testing  backwards  until  the  gradient  again  drops  below  the  threshold 
criteria — ^thus  the  location  of  a  structure  is  bracketed  between  right  and  left  bins.  With 
the  structure  now  bracketed,  the  subroutine  selects  the  mid-point  as  the  location  of  the 
structure  and  stores  the  location  for  further  processing  back  in  the  main  program.  After 
determining  the  location  of  a  structure,  the  subroutine  returns  to  the  main  counter  loop 
and  continues  working  backwards  until  another  qualifying  density  gradient  triggers  the 
bracketing  routine.  At  each  time  step,  the  main  program  compares  the  previous  location 
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of  each  structure  to  the  location  returned  by  the  subroutine.  If  a  structure’s  location  has 
changed,  the  program  updates  the  previous  location  and  writes  the  new  location  to  a  data 
file.  This  procedure  can  track  up  to  eight  shocks  and  discontinuities. 

4.1»3.  Model  Adaptation.  The  MacCormack  algorithm  is  typically  employed  to 
find  numerical  solutions  to  the  Riemann,  or  shock  tube  problem.  In  a  Riemann  problem, 
gases  in  a  cylinder  are  separated  into  two  isolated  regions  by  a  thin,  non-permeable 
membrane  and  maintained  at  different  pressures,  densities,  and  temperatures  as  shown  in 

Figure  9(a).  The  shock  tube  experiment  begins  when  the  membrane  is  ruptured 
suddenly  and  the  two  gases,  initially  at  rest,  are  allowed  to  freely  flow.  The  usual  result 
is  the  formation  of  a  shock  wave  propagating  from  the  region  of  higher  pressure  into  the 
lower  pressure  region,  and  an  expansion  of  lower  pressure  in  the  opposite  direction  as 
depicted  in  Figure  9  (b). 
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(b) 

Figure  9.  Shock  tube  problem. 


Unlike  the  shock  tube  problem,  plasmas  in  the  solar  wind  and  the  magnetosheath 
are  seldom  at  rest.  To  realistically  simulate  the  magnetosheath  system  and  the 
interactions  resulting  from  collision  with  the  bow  shock,  the  boundary  conditions  must  be 
adapted  to  permit  a  steady  flow  into,  and  out  of,  the  model.  The  problem  of  simulating  a 
stationary  bow  shock  at  the  boundary  between  the  solar  wind  and  the  magnetosheath 
requires  a  steady  supersonic  solar  wind  flow  into  the  shock. 

Tascione  compared  the  bow  shock  to  the  aerodynamic  stand-off  shock  that  forms 
in  front  of  a  blunt  obstacle  in  supersonic  wind  tunnel  experiments  (1994:57);  in  fact,  data 
from  Parks  suggests  that  the  bow  shock  maintains  a  fairly  constant  stand-off  distance 
from  the  magnetopause  of  about  two  earth  radii  (1996:501-502).  Figure  10  is  a 
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representation  of  the  earth’s  stationary  bow  shock  characterized  by  a  jump  from  solar 
wind  pressure,  po,  to  magnetosheath  pressure,  pi. 


M, 


Figure  10.  The  bow  shock 

The  bow  shock  will  move  or  remain  stationary  depending  on  what  reference  frame  is 
chosen.  In  Figure  10  two  velocities  are  depicted  on  the  left  side  of  the  bow  shock — ^the 
solar  wind  flow  speed,  Uo,  and  the  bow  shock  propagation  speed,  Vi  =  0,  in  the 
absolute — or  earth-centered — ^reference  frame.  Transforming  to  the  rest  frame  of  the 
solar  wind,  uo  =  0,  and  the  bow  shock  would  move  to  the  left  with  a  shock  speed  Vj  =  Mj 
ao  =  -Uo,  where  Mi  is  the  bow  shock  Mach  number  determined  by  the  upstream  solar 
wind  sound  speed  and  the  flow  speed  into  the  shock.  The  case  in  which  Uo  =  0  and  the 
bow  shock  propagates  to  the  left  is  a  classic  Riemaim  problem  of  a  shock  wave 


56 


propagating  into  a  plasma  at  rest.  Observations  indicate  that  the  bow  shock  is  relatively 
stationary,  i.e.  Vj  =  0,  in  the  absolute  reference  frame,  requiring  that  the  solar  wind  flow 
supersonically  into  the  shock,  and  sub-sonically  behind  the  shock  at  the  magnetopause 
flow  velocity  Ui,  in  order  to  maintain  its  structure  and  location.  Thus,  in  order  to 
maintain  a  stationary  shock  in  a  flowing  plasma,  the  model  was  adapted  to  incorporate 
steady  plasma  flow  into  and  out  of  the  boimdaries. 

A  deficiency  in  the  algorithm  was  the  inability  to  simulate  a  stationary  shock 
wave  in  a  flowing  plasma.  It  was  found  that  a  stationary  shock  would  develop  a 
temperature  trough,  or  a  sharp  decrease  in  temperatures  within  a  few  bins  immediately 
upstream  of  a  shock  front.  See  Figure  1 1 . 
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Figure  11.  Temperature  trough. 


57 


Formation  of  a  temperature  trough  seemed  to  occur  only  for  stationary  shock  waves 
which,  in  the  fluid  reference  frame,  propagate  against  plasma  flow  (see  Figure  10).  The 
magnitude  of  the  drop  in  temperature  at  the  shock  front  also  seemed  to  have  a 
dependence  on  the  magnitude  of  the  stationary  shock,  which  is  determined  by  the  speed 
of  the  upstream  plasma  flowing  into  the  shock.  Temperature  drops  appeared  greater  for 
higher  Mach  numbers:  in  fact,  formation  of  a  temperature  trough  was  not  observed  until 
an  attempt  was  made  to  simulate  a  stationary  Mach  5  shock.  Although  too  few 
simulations  were  made  to  draw  any  definite  conclusions,  it  is  believed  that  formation  of 
temperature  troughs  may  be  a  result  of  numerical  instabilities  in  the  code  arising  from 
shock  waves  propagating  against  the  flow  of  a  plasma.  Figure  1 1  was  produced  from 
data  generated  by  the  MacCormack  algorithm  attempting  to  model  a  stationary  bow 
shock  in  a  plasma  characterized  by  the  solar  wind  parameters  found  in  Table  1,  and  a 
flow  velocity  of  280  km/sec.  The  Mach  number  of  the  bow  shock,  Mj,  was  determined 
to  be  8.44  from  the  definition  of  Mach  number,  M  =  V/ao,  based  on  a  shock  propagation 
speed  of -280  km/sec  and  a  solar  wind  sound  speed  of  33  km/sec. 

The  purpose  of  FCT  is  to  reduce  numerical  oscillations  resulting  from  steep 
parameter  gradients  in  a  finite-difference  approximation  algorithm.  It  was  found  that 
multiplying  the  DF  terms  in  the  FCT  subroutine  by  a  factor  of  1.5;  corrected  temperature 
trough  formation  associated  with  a  stationary  shock  front.  Figure  12  shows  the  effect  of 
multiplying  the  diffusive  terms  of  the  FCT  subroutine  by  a  DF  coefficient  of  1.5,  all  other 
variables  and  inputs  to  the  algorithm  were  exactly  the  same  as  in  the  model  run  that 
produced  Figure  1 1 .  The  simulation  yielded  a  shock  front  that  remained  imchanged  from 
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the  one  depicted  in  Figure  12  through  500  iterations  of  the  algorithm.  Adjusting  the 
amount  of  diffusive  flux  added  to  the  raw  numerical  solution  resulted  in  a  steady, 
stationary  simulation  of  the  bow  shock  in  a  supersonic  plasma  flow. 
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Figure  12.  Flux  corrected  stationary  shock. 

Modeling  the  magnetopause  as  a  stationary  tangential  discontinuity  proved 
impossible  using  a  one  dimensional  code.  The  difficulty  in  modeling  the  magnetopause 
in  a  flowing  plasma,  using  a  one  dimensional  algorithm,  arises  from  the  violation  of  the 
conservation  of  mass  flux  at  the  tangential  discontinuity  where  the  flow  velocity  must  be 
equal  to  zero.  Conservation  of  mass  flux  states  simply  that  the  amount  of  mass  flowing 
into  a  model  bin  must  exactly  equal  the  amoimt  exiting  from  the  opposite  side — what 
comes  in  must  go  out.  In  a  two  or  three  dimensional  model  this  difficulty  does  not  exist: 
if  the  X  component  of  velocity  goes  to  zero,  then  the  y  and  z  velocity  components 
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increase  to  maintain  conservation  of  mass  flux  by  carrying  away  mass  that  would 
otherwise  build  up  in  the  system. 

To  overcome  the  difficulty  of  simulating  interactions  involving  a  stationary 
magnetopause  using  a  one-dimensional  code,  the  problem  was  broken  up  into  two  steps. 
The  first  step  numerically  solved  the  ISW-bow  shock  interaction  separately,  with  the 
magnetopause  removed  fi'om  the  problem,  for  the  velocities  of  the  two  new  shock  waves 

53  and  S4.  Taking  the  magnetopause  out  of  the  model  does  not  affect  the  behavior  S3  or 

54  because  shock  waves  travel  supersonically  and  are  unaffected  by  changes  in  fluid 
conditions  upstream. 

The  second  step  was  to  numerically  solve  the  fluid  equations  for  the  interaction  of 
S4  with  a  tangential  discontinuity  in  a  magnetosheath  where  the  flow  velocity  is  zero.  It 
is  relatively  easy  to  maintain  a  stationary  discontinuity  in  a  fluid  at  rest  because  the  mass, 
momentum,  and  energy  fluxes,  p  u,  p  u^,  and  e  u  respectively,  are  zero  and  are  identically 
conserved.  All  that  is  required  is  to  vary  temperature,  density,  and  magnetic  field  in  such 
a  way  that  total  the  pressure  is  conserved  across  the  discontinuity,  ensuring  that  a 
pressure  gradient  will  not  form  at  the  boundary  and  induce  the  fluid  to  begin  moving. 

The  propagation  of  S4  is  imaffected  when  the  upstream  plasma  is  at  rest  because  the 
Mach  number,  M4,  and  hence  the  propagation  speed  relative  to  the  plasma,  depends  on 
the  shock  strength  or  the  compressibility  (Equations  (12)  or  (1 1)  for  the  gas-dynamic 
case.  Equations  (20)  or  (19)  for  the  MHD  case).  Propagation  speeds  in  the  absolute 
reference  frame  will  differ  by  only  a  constant  from  the  speed  in  the  fluid  reference  frame. 
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4.2.  Model  Validation 


Having  developed  a  numerical  algorithm  for  a  one-dimensional,  explicit,  time 
dependent  solution  of  the  MHD  fluid  conservation  equations,  the  algorithm  was  tested 
against  known  analytical  solutions  to  see  if  it  would  behave  as  expected. 

The  first  test  was  the  most  straightforward:  a  simple  gas-dynamic  Mach  3.5  shock 
propagating  through  a  plasma  at  rest.  In  this  test  the  unshocked  upstream  plasma  was 
characterized  by  the  ambient  solar  wind  parameters  for  pressure,  density,  and 
temperature  listed  in  Table  1.  The  analytic  solution  for  the  normalized  shock 
propagation  speed  in  terms  of  Mach  number  and  sound  speed  is. 


where, 

V  =  shock  propagation  speed, 

M  =  shock  Mach  number, 
ao  =  solar  wind  sound  speed, 

Vth  =  solar  wind  thermal  velocity  =  ( k  To/m)'^. 

The  analytic  solution  is  normalized  to  the  solar  wind  thermal  velocity,  Vo,,  to  be 
consistent  with  the  normalization  of  model  parameters  (see  parameter  transformations, 
page  47).  For  a  Mach  3.5  shock  propagating  in  a  stationary  solar  wind  having  a  sound 
speed  of  33  km/sec  and  normalized  to  a  thermal  velocity  of  1 8  km/sec  (from  Table  1),  Eq 
(73)  yielded  a  shock  speed  of  6.3901,  shown  as  the  horizontal  line  in  Figure  13. 
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Positions  for  the  same  Mach  3.5  shock  wave  were  determined  numerically  by  the 
algorithm  and  were  analyzed  using  a  moving  average  least  square  fit  to  the  data  taken 
over  a  twenty-one  point  interval — ^ten  data  points  on  either  side  of  the  point  at  which  the 
velocity  was  to  be  found — in  order  to  estimate  shock  propagation  speeds  from  the  inverse 
of  the  linear  coefficient  of  the  fitted  line.  The  oscillatory  line  in  Figure  9  is  a  plot  of 
least-square  averaged  velocities  for  the  Mach  3.5  shock.  The  mean  averaged  shock 
velocity  was  6.3914,  agreeing  to  within  less  than  a  tenth  of  a  percent  of  the  analytic 
solution. 


Figure  13.  Mach  3.5  shock  velocity 

The  next  test  was  designed  to  validate  the  ability  of  the  numerical  model  to 
simulate  a  shock  wave-discontinuity  interaction.  In  this  case,  a  Mach  3.5  shock 
propagating  in  a  stationary  gas-dynamic  solar  wind  impinges  on  a  thermal  discontinuity. 
Temperature  in  the  “hot”  region  was  increased  by  a  factor  of  ten  over  the  solar  wind 
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temperature.  To  ensure  continuity  of  pressure  across  the  discontinuity,  the  density  was 
correspondingly  decreased  to  one  tenth  of  the  solar  wind  density.  Applying  the  shock 
jump  equations  (Equations  (1 1)  -  (14))  to  a  Mach  3.5  shock  wave  propagating  into  an 
unshocked  solar  wind  determined  the  plasma  parameters  behind  the  shock.  The  plasma 
parameters  for  each  test  region,  normalized  to  ambient  solar  wind  parameters,  are 
presented  in  Table  4. 


Table  4.  Gas-dynamic  shock-discontinuity  validation  palsma  parameters. 


Parameter 

Shocked 

Region 

Stationary 
Solar  Wind 

Thermalized 

Region 

P 

15.06 

1.00 

1.00 

P 

3.21 

1.00 

0.10 

T 

4.69 

1.00 

10.0 

u 

4.40 

0.00 

0.00 

a 

3.95 

1.83 

5.77 

Recall  that  a  shock  impinging  on  a  discontinuity  results  in  a  new  transmitted 
shock  wave,  a  reflected  rarefaction  wave,  and  a  new  contact  discontinuity  as  depicted  in 
Figure  5.  The  Mach  number  of  the  transmitted  shock.  Me,  was  determined  analytically 
by  first  solving  the  following  equation  for  yi. 


0-^)(yi-i) _ ^ 

./O  +  pKy^Tp)  Y-1 


Vy.; 


-1 


yi 


i+ii+y.y'' 

ii+yi  ) 


jW  0-^)(y2-i) 

V  p2  Vo+p)(y2+p) 


(74) 
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with. 


^  =  (Y-1)/(Y+1), 

y2  =  shock  strength  in  the  thermalized  region, 
yi  s  shock  strength  of  the  shock  in  the  solar  wind. 

Pi  =  density  in  the  solar  wind, 

P2  =  density  in  the  thermalized  region  (Wright,  1961:79); 
and  then  by  setting  y2  equal  to  the  gas-dynamic  shock  strength  (Eq  (12)).  With 
known,  the  speed  of  the  transmitted  shock.  Ye,  was  then  found  from  Eq  (73),  and  the 
flow  speed  behind  the  transmitted  shock,  14,  was  easily  solved  using  Eq  (13).  Flow 
velocity  behind  the  transmitted  shock  is  of  interest  because  it  is  the  speed  at  which  the 
contact  surface  will  move. 

The  numerical  algorithm  was  executed  for  a  Mach  3.5  shock  propagating  in  a 
stationary  solar  wind  plasma  in  which  the  same  thermal  discontinuity — i.e.  T  =  10.0  To 
and  p  =  0. 1  po  —was  initialized.  The  model  tracked  the  positions  of  the  transmitted 
shock  and  the  contact  discontinuity,  but  failed  to  locate  the  rarefaction  wave,  except 
perhaps  at  the  onset  of  the  modeled  interaction,  as  illustrated  by  the  x-t  plot  of  shock 
position  data.  The  locator  subroutine  failed  to  detect  the  rarefaction  wave  because 
changes  in  density  are  spread  out  over  many  bins,  thus  decreasing  density  gradients 
within  the  wave  below  the  locator  threshold  criteria. 
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Figure  14.  Shock-discontinuity  x-t  diagram. 


Analytic  and  numerical  mean  averaged  velocities  for  and  ug  are  shown  in  Table  5. 
Again  the  numerical  model  was  within  less  than  a  percent  of  the  analytical  solution. 


Table  5.  Gas-dynamic  shock-discontinuity  validation  results. 


Speed 

Analytic 

Numerical 

% 

Results 

Result 

Error 

V6 

11.82 

11.84 

-0.17 

U6 

6.752 

6.800 

-0.71 

The  last  test  of  the  model  using  the  gas-dynamic  governing  fluid  equations  is  a 
shock-shock  interaction  representing  the  ISW-bow  shock  collision.  Here  a  Mach  3.5 
interplanetary  shock  propagating  through  a  solar  wind  flowing  at  280  km/sec  impinges  on 
a  stationary  bow  shock  resulting  in  two  shock  waves,  S3  and  S4,  with  a  contact 


discontinuity,  Cj,  between  them.  Initial  nonnalized  plasma  parameters  are  given  in 
Table  6. 


Table  6  .  Gas-dynamic  shock-shock  validation  plasma  parameters. 


Parameter 

Shocked 

Region 

Solar  Wind 

Magnetosheath 

Region 

P 

15.06 

1.00 

88.70 

P 

3.21 

1.00 

3.84 

T 

4.69 

1.00 

23.11 

u 

19.80 

15.40 

4.01 

a 

3.95 

1.83 

8.78 

Grib’s  equations  ( Equations  (39)  and  (40))  describing  the  continuity  of  flow 
velocity  and  pressure  across  Cj  are  solved  simultaneously  to  obtain  values  for  M3  and 
M4;  3.60  and  1.54  respectively.  The  propagation  speeds;  of  S3  and  S4,  V3  and  V4,  are 
found  after  substituting  M3  and  M,  into  Eq  (73),  and  Eq  (35)  with  M3  results  in  a 
determination  of  the  speed  U3,  and  thus  the  speed  of  Ci. 

The  analytic  and  numerical  solutions  for  the  velocities  of  the  shock  waves  and  the 
contact  discontinuity  are  presented  in  Table  7. 


Table  7.  Gas-dynamic  shock  shock  validation  results. 


Speed 

Analytic 

Result 

Numerical 

Result 

% 

Error 

V3 

5.55 

5.60 

0.90 

V4 

17.57 

17.86 

1.65 

U3 

9.92 

9.91 

0.10 
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After  demonstrating  the  ability  of  the  algorithm  to  produce  reasonable 
approximations  of  shock  interactions  consistent  with  analytic  solutions  for  a  gas-dynamic 
plasma,  the  next  step  in  the  validation  procedure  was  to  repeat  the  tests  to  show  that  the 
model  also  reasonably  approximated  shock  interactions  and  movement  for  a  MHD 
plasma.  As  in  the  validation  of  the  gas-dynamic  model,  the  first  test  of  the  MHD 
algorithm  was  to  compare  the  numerical  solution  for  a  Mach  3.5  shock  wave  propagating 
through  a  stationary  solar  wind  plasma  and  compare  the  result  to  the  expected  analytical 
result,  V=6.3900,  obtained  from  Eq  (73).  Analysis  of  the  numerical  data  yields  a  mean 
averaged  velocity  of  6.39104,  within  0.016%  of  the  analytic  value. 

Following  the  same  procedure  that  was  used  to  validate  the  gas-dynamic 
algorithm,  the  next  test  of  the  MHD  model  was  an  attempt  to  simulate  a  shock- 
discontinuity  interaction  within  a  stationary  solar  wind  plasma.  The  attempt  to  model  the 
Mach  3.5  shock-discontinuity  interaction  failed.  An  attempt  was  made  to  substitute  a 
Mach  1.5  shock  in  the  problem  with  similar  results.  A  last  attempt  to  achieve  positive 
results  involved  switching  to  a  Mach  8.0  shock  propagating  through  a  magnetosheath 
plasma  in  pressure  balance  with  a  more  realistic  magnetopause.  Plasma  parameters  for 
the  magnetosheath  and  magnetosphere  used  in  this  simulation  are  listed  in  Table  8. 

Table  8.  MHD  magnetosheath  and  magnetopause  parameters. 


Parameter 

Region  1 

Magnetopause 

Density ,  p  (x  10'^°  kg/m^) 

6.85 

0.68 

Temperature,  T  (x  10“*  °K) 

88.6 

8.86 

Flow  Velocity,  u  (km/sec) 

75.0 

0.00 

Magnetic  Field,  B  (nT) 

13.1 

51.6 

Pressure,  P  (x  10‘"  Pa) 

1.08 

0.01 

Sound  Speed,  a  (km/sec) 

156 

49.4 

Plasma  Beta 

14.8 

0.009 
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Parameters  in  Table  8  are  chosen  to  be  consistent  with  Grib’s  treatment  of  the 
magnetopause. 

The  Mach  8.0  simulation  succeeded  in  resolving  a  transmitted  shock,  a  contact 
discontinuity,  and  suggestions  of  a  rarefaction  wave  from  location  data,  but  velocity 
approximations  from  the  model  differed  significantly  (~  14-20%)  from  analytic  solutions. 
The  result  of  the  Mach  8.0  simulation  are  given  in  Table  9  along  with  analytic  solutions 
and  solutions  from  Grib  (1979:5912)  of  a  similar  calculation.  The  analytic  solutions 
were  derived  from  the  MHD  continuity  of  plasma  flow  velocity  and  total  pressure 
(Equations  (45)  and  (46))  for  a  stationary  solar  wind  plasma  with  a  magnetic  field  of 
3.5  nT.  Grib  applies  the  same  continuity  equations  to  a  solar  wind  plasma  in  which  the 
magnetic  field  is  2.5  nT  (1979:5912). 

Table  9.  MHD  shock-discontinuity  validation  results. 


Speed 

Analytical 

Numerical 

Grib’s 

Solution 

Results 

Solution 

V6 

16.04 

18.35 

15.69 

U6 

8.88 

10.66 

12.55 

The  failure  of  the  model  to  simulate  the  shock-discontinuity  interaction  for  the 
MHD  plasma  is  disturbing.  However,  failure  in  this  one  particular  simulation  should  not 
be  taken  as  evidence  that  the  model  does  not  reasonably  approximate  these  kinds  of 
interactions,  especially  in  light  of  other  successful  simulations  of  shock  wave 
propagation  and  shock-shock  interactions.  It  is,  rather,  an  indication  that  further  work  is 
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needed  to  refine  the  representation  of  a  stationary  discontinuity  and  the  equations  that 
determine  the  analytic  solutions. 


On  a  positive  note,  the  last  simulation  to  test  the  ability  of  the  model  to  simulate 
shock-shock  interactions  in  a  MHD  fluid  was  successful.  A  Mach  3.5  shock  was 
propagated  into  a  non-stationary  solar  wind  plasma  flowing  at  280  km/sec  and  allowed  to 
collide  with  a  stationary  bow  shock.  The  model  handled  this  particular  simulation  very 
well.  The  location  of  each  of  the  resultant  shock  wave  pair,  and  the  contact  discontinuity 
between  them,  was  tracked  quite  effectively  by  the  program  as  seen  in  the  figure  below. 


Figure  15.  Numerical  shock-shock  interaction. 
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The  analytic  solutions  for  shock  speed  were  derived  as  before:  the  equations  of 
continuity  of  flow  velocity  and  total  pressure  (Equations  (45)  and  (46))  across  the  contact 
discontinuity  Ci  were  solved  simultaneously  to  determine  M3  and  M4.  The  upstream 
flow  into  S3  is  a  plasma  previously  shocked  by  the  initial  Mach  3.5  shock  wave;  the 
upstream  sound  speed,  ^2,  is  found  by  applying  the  shock  jump  equations  (Equations  (19) 
-  (23))  with  M  ==  3.5  to  solar  wind  parameters  and  then  substituting  the  resultant  pressure 
and  density  into  equation  (2).  Similarly,  the  sound  speed  upstream  from  S4  is  just  the 
magnetosheath  sound  speed  and  Eq  (73),  applied  to  the  two  new  shock  waves  results  in 
the  normalized  propagation  speeds.  The  speed  of  Ci  is  the  flow  speed  of  the  plasma 
between  the  two  shocks  determined  by  either  the  left  or  right  sides  of  Eq  (43).  The 
analytical  solutions  are  shown  with  the  numerical  results  in  Table  10. 

Table  10.  MHD  shock-shock  validation  results. 


Speed 

Analytic 

Results 

Numeric 

Results 

% 

Error 

V3 

5.15 

4.89 

5.07 

V4 

17.32 

17.58 

1.50 

U3 

9.45 

9.59 

1.48 

Comparing  the  results  presented  in  Tables  7  and  10  for  the  gas-dynamic  and 
MHD  solutions  of  the  shock-shock  interactions  it  is  seen  that  magnetic  effects  within  the 
model  slowed  V4  and  V3  slightly  (by  ~  two  percent  and  twelve  percent  respectively,  in 
this  case),  and  slightly  increases  the  velocity  of  the  contact  discontinuity.  This 
observation  is  consistent  with  Shen’s  estimate  that  the  magnetic  field  slowed  the 
propagation  speed  of  S4  and  S3  by  about  ten  percent  (Shen  and  Dryer,  1972:4635). 
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4.3  Model  Initialization  and  Execution 


The  final  section  in  this  chapter  is  a  discussion  of  the  methodology  of  model 
initialization ,  execution,  and  data  processing.  Plasma  pressure,  density,  temperature, 
magnetic  field  (for  MHD  fluids  only),  and  flow  velocity  must  be  initialized  by  the  user  at 
each  bin  in  the  numerical  model.  In  addition,  each  parameter  must  be  normalized  to  an 
ambient,  unshocked  plasma  determined  from  number  density,  temperature,  magnetic 
field,  and  flow  velocity  supplied  by  the  user  as  input  to  the  model.  The  program  assists 
in  the  initialization  process  by  determining  appropriate  normalized  parameters  for  the 
bow  shock  and  ISW  based  upon  the  shock  jump  equations  and  a  Mach  number  for  the 
ISW,  which  is  also  provided  by  the  user,  but  parameter  values  at  each  bin  must  be  set 
manually  from  inside  the  code.  Other  input  variables  used  by  the  algorithm  include  the 
scale  length  of  the  physical  problem  and  the  number  of  bins  desired,  which  determines 
the  spatial  size  of  the  grid.  The  user  selects  a  termination  criteria:  a  specified  number  of 
iterations,  or  at  a  specified  end  time.  Initial  locations  in  normalized  x  coordinates  of  the 
ISW,  bow  shock  and  magnetopause  are  also  input  by  the  user. 

The  input  parameters  no.  To,  Bo,  and  M2  used  to  validate  the  model  were  selected 
to  coincide  with  Grib,  et.  al.  In  their  paper,  they  chose  three  values  for  the  incident  shock 
wave  Mach  number,  M2  =  1.5, 3.5,  and  8.0,  as  a  representative  range  from  low  to  high 
values. 

Once  the  model  was  initialized,  execution  proceeded  fairly  quickly  on  a  desk  top 
computer.  Typical  execution  times  were  about  5  -  20  minutes  depending  on  the  size  of 
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the  grid  and  the  speed  of  the  shock  wave  simulated.  Large  grid  sizes  produced  large  data 
files  and  a  disproportionate  time  was  spent  writing  to  output.  A  gut  feel  choice  was 
made  to  adopt  a  2401  (it  must  be  an  odd  number)  bin  grid  as  a  compromise  between  fine 
spatial  resolution  and  long  simulation  times,  and  a  coarse  grid  with  a  fast  execution  time 
that  might  not  spatially  resolve  shock  fronts  or  density  discontinuities.  Strong  shocks 
with  a  high  Mach  number,  >5,  propagate  quickly  across  a  model  bin,  requiring  a  very 
short  time  step  in  order  to  prevent  flow  across  more  than  one  model  bin.  A  shorter  step 
time  increased  both  computation  time  and  data  storage  requirements. 

The  ability  to  process  data  files  generated  by  the  model  into  x-t  diagrams^  and 
into  pressure,  density,  temperature,  and  velocity  profiles  that  could  then  be  animated  to 
show  the  time  evolution  of  the  numerical  simulation  was  a  useful  tool  to  understand  both 
the  physical  and  numerical  processes. 

The  development  and  validation  of  an  algorithm  to  numerically  solve  the  one¬ 
dimensional  gas-dynamic  and  MHD  fluid  equations,  and  extraction  of  shock  velocities 
and  locations  from  the  output  data,  has  demonstrated  the  applicability  and  utility  of  this 
analysis  tool  to  the  investigation  of  bow  shock  and  magnetopause  displacement  resulting 
from  an  impinging  interplanetary  shock  wave.  In  the  next  chapter  the  results  of  the 
application  of  the  algorithm  developed  here  will  be  presented  and  discussed. 
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V.  Results 


This  chapter  is  a  presentation  of  the  results  of  a  systematic  application  of  the 
numerical  algorithm  developed  in  this  thesis  to  the  problem  of  determining  initial  bow 
shock  and  magnetopause  displacement  following  interaction  with  an  interplanetary  shock 
wave.  For  both  the  gas-dynamic  and  the  MHD  cases,  three  incident  interplanetary 
shocks— having  Mach  numbers  M2  =  1.5, 3.5,  and  8.0— were  simulated  interacting  with 
the  earth’s  bow  shock.  The  algorithm  was  also  applied  to  the  interaction  of  the 
magnetopause  with  three  different  shock  waves,  S4,  resulting  from  the  ISW-bow  shock 
collision.  Numerical  solutions  for  shock  and  discontinuity  velocities,  derived  from  data 
resulting  from  these  simulations  are  presented  for  comparison  with  normalized  analytic 
velocities  determined  from  the  application  of  the  appropriate  forms  of  the  equations  of 
continuity  across  contact  discontinuities  and  Eq  (73). 

The  first  set  of  results  are  from  the  gas-dynamic  and  MHD  simulations  of  the 
ISW-bow  shock  collision.  The  ambient  plasma  in  the  these  simulations  was  the  quiet 
solar  wind  flowing  at  280  km/sec,  in  which  a  stationary  bow  shock  was  maintained.  The 
normalized  results  are  presented  in  tabular  form  in  Tables  1 1  and  12.  A  small  xx  in  the 
tables  indicates  missing  data  due  to  the  failure  of  the  locator  subroutine  to  detect  small 
gradients  in  the  output  data.  The  symbol  A%  represents  the  percentage  error  between  the 
“true”  analytic  and  the  time-averaged  numerical  solutions,  defined  by. 
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A%  = 


Analytic  -  Numeric 
Analytic 


xlOO. 


Table  1 1.  Gas-dynamic  shock-shock  results. 


Incident 

Mach 

Number 

Type 

Solution 

V3 

V4 

U3 

M2  =  1.5 

Analytic 

2.32 

14.30 

6.14 

Numerical 

2.29 

14.32 

6.05 

A% 

1.29 

-0.14 

1.47 

M2  =  3.5 

Analytic 

5.55 

17.57 

9.94 

Numerical 

5.60 

17.86 

9.91 

A% 

-0.90 

-1.65 

0.30 

II 

00 

o 

Analytic 

7.51 

22.73 

14.97 

Numeric 

7.49 

23.96 

14.51 

A% 

0.26 

-5.41 

-3.10 

Table  12.  MHD  shock-shock  results. 


Incident 

Mach 

Number 

Type 

Solution 

V3 

V4 

U3 

M2  =  1.5 

Analytic 

2.24 

14.30 

6.06 

Numerical 

-1.4 

11.45 

2.01 

A% 

162 

20 

66.8 

M2  =  3.5 

Analytic 

5.15 

17.34 

9.89 

Numerical 

4.89 

17.57 

9.58 

A% 

5.04 

-1.32 

3.13 

P 

00 

II 

Analytic 

7.50 

22.75 

14.83 

Numeric 

7.31 

22.88 

14.64 

A% 

2.53 

-0.57 

1.28 
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Tables  1 1  and  12,  show  that  the  algorithm  closely  approximated  the  velocities 
determined  by  the  anal)^ic  equations  for  V3  and  V4  to  within  about  five  percent  in  all  but 
one  simulation.  An  analysis  of  shock  positions  from  the  MHD  simulation  of  the  M2  = 

1 .5  shock  indicated  a  negative  value  for  V3  and  a  lower  value  of  V4.  Movement  of  S3  in 
the  negative  (left)  direction  indicates  that  the  propagation  speed,  V3=  M3  a2,  of  the  new 
bow  shock  into  the  upstream  ISW  plasma  is  greater  than  the  flow  velocity  U2  since,  in  the 
rest  reference  frame,  shock  propagation  speed  is  determined  by  U2  -  V3;  this  problem 
could  arise  if  the  solar  wind  flow  velocity,  Uo,  was  mistakenly  set  to  zero  in  the  numerical 
model.  A  stationary  solar  wind  would  also  explain  propagation  speeds  for  Cj  and  V4  less 
than  analytical  predictions  because  the  magnetosheath  flow  velocity  depends  on  the  solar 
wind  flow  speed  from  the  shock  jump  equations — ^zero  magnetosheath  flow  speed 
translates  into  slower  than  expected  V4  and  U3. 

The  second  set  of  data  presented  are  the  result  of  gas-dynamic  simulations  of  the 
S4-magnetopause  interaction.  Because  S4  results  from  the  interaction  of  the  ISW  with  the 
bow  shock,  to  properly  initialize  the  model  required  prior  determination  of  M4,  and  all 
plasma  parameters  associated  with  the  shock,  using  the  equations  of  continuity  for  flow 
velocity  and  pressure  and  the  shock  jiunp  equations.  In  this  simulation,  magnetospheric 
plasma  density  was  assumed  to  be  one-tenth  the  density  in  the  magnetosheath,  consistent 
with  Grib’s  analysis;  to  maintain  a  balance  of  pressure  across  the  tangential  discontinuity, 
a  magnetospheric  temperature  equal  to  ten  times  the  magnetosheath  temperature  was 
required. 
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The  Mach  number  for  the  transmitted  shock.  Me,  was  determined  from  Grib’s 
equations  for  continuity  of  flow  velocity  and  pressure  (Equations  (39)  and  (40))  across 
the  new  magnetopause,  Because  the  plasma  was  assumed  to  be  an  adiabatic  ideal 
gas,  as  and  ps,  the  sound  speed  and  pressure  in  the  rarefaction  wave,  had  the  forms  of 
Equations  (41)  and  (42).  M^  and  a^,  the  magnetospheric  sound  speed  determined  by  the 
initialization  of  the  magnetosphere,  determined  the  shock  speed  when  substituted  into  Eq 
(73). 


Table  13.  Gas-dynamic  shock-discontinuity  results. 


Incident 

Mach 

Number 

Type 

Solution 

V6 

U6 

M2  =  1.5 

Analytic 

30.00 

3.18 

Numerical 

29.96 

3.19 

A% 

0.007 

-0.063 

M2  =  3.5 

Analytic 

34.37 

8.97 

Numerical 

34.43 

8.98 

A% 

-0.19 

-0.15 

0 

00 

II 

Analytic 

41.01 

16.67 

Numeric 

41.05 

16.70 

A% 

-0.11 

-0.20 

The  last  set  of  results  presented  are  the  from  the  M2  =  8.0  MHD  simulations  of 
the  S4-magnetopause  interaction  shown  in  Table  14.  Again,  the  large  discrepancies 
between  analytical  and  numerical  solutions,  and  the  lack  of  usable  data  for  the  M2  =  1.5 
and  3.5  simulations  suggest  a  flaw  in  the  numerical  code  to  model  the  magnetopause  in  a 
MHD  plasma,  errors  in  the  analytic  solutions,  or  a  combination  of  both. 
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Table  14.  MHD  Mach  8.0  shock-discontinuity  results. 


M  =  8.0 
Results 

Type 

Solution 

V6 

U6 

Analytic 

Analytic 

16.04 

8.88 

Numerical 

18.35 

10.66 

A% 

-14.4 

-20.0 

Grib 

Grib 

15.69 

12.55 

Numerical 

18.35 

10.66 

A% 

-17.0 

15.1 

The  good  agreement  between  analytic  and  numerical  solutions  for  velocities  of 
the  two  new  shock  waves  resulting  from  collision  of  an  interplanetary  shock  wave  and 
the  earth’s  bow  shock,  as  well  as  the  gas-dynamic  solutions  to  the  shock-discontinuity 
interaction,  suggests  that  the  algorithm  developed  in  this  thesis  is  sound.  Though 
discouraging,  the  failure  of  the  algorithm  to  closely  approximate  the  analytical  results  for 
the  MHD  shock-magnetopause  simulations  indicate  a  need  for  refinement  of  the  analytic 
solutions  to  eliminate  possible  errors,  and  of  the  initialization  of  plasma  parameters  put 
into  the  model  to  represent  the  magnetopause.  These  results  show  that  the  algorithm 
produces  reasonable  approximations  to  plasma  behavior  consistent  with  an  analytic 
analysis,  primarily  by  Grib,  et.  al.,  which  describe  the  physical  processes  modeled  in 
these  simulations:  a  logical  next  step  is  a  thorough  validation  of  the  model  with  respect  to 
real  observations  of  the  bow  shock  and  magnetopause  motions.  The  kind  of  validation 
required  is  becoming  more  feasible  as  the  number  of  space  platforms  and  observation 
techmques  continues  to  expand  the  database  of  bow  shock  and  magnetopause  behaviors. 
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VI.  Conclusions  and  Recommendations 


This  thesis  had  four  objectives.  The  first  was  to  develop  an  understanding  of  the 
physical  processes  associated  with  the  displacement  of  the  earth’s  bow  shock  and 
magnetopause  initiated  by  collision  of  an  ISW  with  the  magnetosheath  system.  The 
second  objective  was  the  development  of  a  one-dimensional  numerical  simulation  of  the 
ISW-magnetosheath  interactions  based  upon  gas-dynamic  and  MHD  fluid  equations  to 
determine  the  displacements  and  velocities  of  the  bow  shock  and  the  magnetopause.  A 
third  objective  was  to  review  and  understand  prior  analytic  treatments  of  the  ISW- 
magnetosheath  interactions,  providing  a  basis  for  comparison  with  numerical  results. 
And  the  fourth  objective  was  to  evaluate  the  performance  of  the  numerical  model  by 
comparing  simulation  results  with  analytic  solutions.  In  terms  of  development  and 
validation  of  an  algorithm  to  solve  the  one-dimensional  fluid  equations,  and  the 
application  of  the  algorithm  to  simulations  of  the  ISW-magnetosheath  interactions,  this 
thesis  has  accomplished  each  of  its  objectives. 

Based  on  the  validation  of  the  algorithm  and  simulations  of  the  ISW-bow  shock 
interactions,  the  development  of  the  governing  fluid  equations  and  their  incorporation 
into  a  numerical  algorithm  to  model  these  interactions  is  sound.  A  problem  simulating 
the  shock-magnetopause  interactions  does  not  cast  doubts  on  the  basic  soundness  of  the 
model;  it  does  indicate  that  the  MHD  treatment  of  the  shock-magnetopause  interaction 
requires  further  development. 
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During  the  course  of  this  thesis  project,  several  ideas  aimed  to  improve  the  model 
and  enhance  its  utility  in  the  simulation  of  bow  shock  and  magnetopause  interaction  were 
developed  and  are  presented  as  recommendations  to  those  interested  or  who  may  wish  to 
continue  this  work. 

The  first  recommendation  is  to  expand  the  one  dimensional  model  into  a  two  or 
three  dimensional  code  in  order  to  better  simulate  the  entire  interaction  of  an  ISW  with 
the  magnetosheath  system.  A  multidimensional  code  would  permit  a  seamless 
simulation  without  the  necessity  of  dividing  the  problem  into  two  separate  steps.  A 
multidimensional  code  would  also  allow  a  more  realistic  simulation  of  a  magnetosheath 
in  which  the  flow  velocity  at  the  stagnation  point  equaled  zero  without  violating  the 
conservation  of  mass,  momentum,  and  energy  fluxes.  A  seamless  simulation  of  the  ISW- 
magnetosheath  interactions  would  conceivably  permit  numerical  simulation  of  multiple 
rarefaction  wave  reflections  between  the  advancing  bow  shock  and  magnetopause  which 
are  proposed  by  Grib,  et.  al.  as  the  mechanisms  that  eventually  slow  and  reverse  the 
magnetopause  during  magnetospheric  compression. 

A  second  recommendation  is  that  the  initialization  of  the  magnetosphere  in  the 
model  be  verified  and  corrected  if  necessary,  and  that  the  simulations  of  the  shock- 
magnetopause  interactions  be  re-accomplished. 

A  third,  and  final,  recommendation  is  that  this  algorithm  be  validated,  as  the  data 
becomes  available,  by  comparison  of  numerical  solutions  to  observations  of  actual  bow 
shock  and  magnetopause  displacements.  The  validation  may  suggest  ways  that  the 
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algorithm  can  be  improved  or  applied  to  other  problems  of  interest  in  plasma  or  shock 
physics. 


80 


Vn.  Bibliography 


Courant,  R.  and  K.O.  Friedrichs.  Supersonic  Flow  and  Shock  Waves.  New  York: 
Interscience  Publishers,  1948. 

Fletcher,  C.A.J.  Computational  Techniques  for  Fluid  Dynamics  Vol.  II:  Specific 
Techniques  for  Different  Flow  Categories.  2d  ed.  Berlin:  Springer- Verlag,  1991. 

Grib,  S.A.,  et.  al.  “Interaction  of  Interplanetary  Shock  Waves  with  the  Bow  Shock- 
Magnetopause  System,”  Journal  of  Geophysical  Research.  84.  No.  Al  0:5907-5921 
(October,  1979). 

Hilbun,  W.  A  One  Dimensional  MacCormack  Algorithm  with  FCT.  Unpublished 
software.  Air  Force  Institute  of  Technology,  Da34on  (1997). 

Hudson,  M.K.,  et.  al.  “Simulation  of  Proton  Radiation  Belt  Formation”,  Geophysical 
Research  Letters.  20  NO.  3:291-294  (February,  1995). 

Kivelson,  M.  and  C.  Russell.  Introduction  to  Space  Physics.  New  York:  Cambridge 
1995. 

Landau,  L.  D.  andE.  M.  Lifshitz.  Fluid  Mechanics.  London:  Pergamon  Press,  1959. 

Li,  Xinlin,  et.  al.  “Simulation  of  the  Prompt  Energization  and  Transport  of  Radiation 
Belt  Particles  During  the  March  24, 1991  SSC,”  Geophysical  Research  Letters.  20. 

No.  22:2423-2426  (November,  1993). 

Parks,  George  K.  Physics  of  Space  Plasmas:  An  Introduction.  New  York;  Addison- 
Wesley,  1996. 

“SOHO  Orbit  Schematic.”  Excerpt  from  the  SOHO  Home  Page,  n.  page.  WWWeb 
http://sohowww.nascom.nasa.gov/gifrhalo_orbit.gif  6  November,  1997. 

Tascione,  T.  F.  Introduction  to  the  Space  Environment.  2d  ed.  Malabar:  Krieger 
Publishing,  1994. 

Toth,  G.  and  D.  Odstr6il.  “Comparison  of  Some  Flux  Corrected  Transport  and  Total 
Variation  Diminishing  Numerical  Schemes  for  Hydrodynamic  and 
Magnetohydrodynamic  Problems,”  Journal  of  Computational  Physics.  128.:82-1Q0 
(1996). 

Shen,  W.W.  “Interaction  of  Interplanetaiy  MHD  Shock  Waves  with  the  Magnetopause,” 
Astrophysics  and  Space  Science.  24:51-64  /Febmarv  1973). 


81 


Shen,  W.W.  and  M.  Dryer.  “Magnetohydrodynamic  Theory  for  the  Interaction  of  an 
Interplanetary  Double-Shock  Ensemble  with  the  Earth’s  Bow  Shock,”  Journal  of 
Geophysical  Research.  77,  No.  25:4627-4644  (September,  1972). 

Vvedensky,  D.  D.  Partial  Differential  Equations  With  Mathematica.  London:  Addison- 
Weseley,  1993. 

Wright,  J.K.  Shock  Tubes.  London:  Methuen,  1961. 


82 


Vita 


Capt  William  A.  Olson  was  bom  on  13  October  1963  in  Tooele,  Utah.  He 
graduated  from  Tooele  High  School  in  1982  and  entered  undergraduate  studies  at 
Brigham  Young  University  in  Provo,  Utah  in  August  1984.  Upon  graduation  in  1988 
with  a  Bachelor  of  Science  degree  in  Physics,  he  entered  Officer  Training  School  in  San 
Antonio,  Texas  and  was  commissioned  a  2nd  Lieutenant  in  May  1989.  For  the  next  year 
he  attended  the  University  of  Oklahoma  in  Norman,  Oklahoma  where  he  completed  the 
Basic  Meteorology  Program  in  May  1990. 

Capt  Olson’s  first  assignment  began  in  June  1990  attached  to  the  XVIII  Airborne 
Corps  at  Fort  Bragg,  North  Carolina  as  a  Staff  Weather  Officer  (SWO).  After  completing 
parachute  training,  Capt  Olson  was  sent  to  Saudi  Arabia  as  part  of  Operations  Desert 
Shield  and  Desert  Storm  where  he  served  with  the  XVIII  Airborne  Corps  and  the  24th 
Infantry  Division  (Mechanized)  in  northern  Saudi  Arabia  and  southern  Iraq.  Following 
the  Gulf  War  he  returned  to  Fort  Bragg  to  serve  as  XVIII  Airborne  Corps  and  82nd 
Airborne  Division  SWOs  until  1993. 

In  February  1993  Capt  Olson  was  reassigned  to  US  Central  Command  at  MacDill 
AFB,  Florida.  While  at  MacDill  AFB,  he  earned  a  Master  of  Aeronautical  Science 
degree  in  Management  from  Embry-Riddle  Aeronautical  University.  In  May  1996  he 
entered  the  School  of  Engineering  Physics,  Air  Force  Institute  of  Technology. 

Permanent  Address:  606  N  Broadway  Ave 
Tooele,  UT  84074 


83 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704^0188 


FHibk  reporting  burden  for  this  collwtion  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing 
tw  collection  of  mtormation  Sand  cornmerits  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  Information 
Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  222024302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503 


1.  AGENCY  USE  ONLY  2.  REPORT  DATE  3.  REPORT  TYPE  A1 

_ DEC  97 _ 

~4.  TITLE  AND  SUBTITLE  - 

DISPLACEMENT  OF  THE  EARTH’S  BOW  SHOCK  AND  MAGNETOPAUSE 
DUE  TO  AN  IMPINGING  INTERPLANETARY  SHOCK  WAVE 


3.  REPORT  TYPE  AND  DATES  COVERED 

_ MASTER’S  THESIS 

I  5.  FUNDING  NUMBERS 


6.  AUTHOR(S) 


William  A.  Olson,  Captain,  USAF 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Institute  of  Technology 
2750  P  Street 

Wright-Patterson  AFB,  OH  45433-7765 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESSjES) 

PL/GPS 

Hanscom  AFB,  MA  01731-5000 

11.  SUPPLEMENTARY  NOTES 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFIT/GAP/ENP/97D-08 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


,  12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  unlimited 

13  ABSTRACT  (Maximum  280  words)  “““  ~  - - - - 

Interplanetary  shock  waves  (ISWs)  propagating  through  the  solar  wind  can  "collide"  with  the  earth’s  bow  shock,  resulting 
in  a  series  of  new  shocks,  contact  discontinuities,  and  rarefaction  waves  which  interact  to  effectively  move  the  bow  shock 
and  magnetopause  toward  the  earth.  A  one-dimensional  MacCormack  predictor-corrector  algorithm  with  Flux  Corrected 
Transport  (FCT)  was  developed  to  model  the  ISW-bow  shock  and  magnetopause  interactions,  and  to  mumerically  predict 
their  propagation  speeds  after  collision.  Analytic  relationships  for  the  Mach  numbers  and  propagation  speeds  of  the 
generated  shock  waves  and  contact  discontinuities  were  used  to  validate  the  model  and  to  compare  numerical  results.  In  both 
the  gas-dynamic  and  magnetohydrodynamic  (MHD)  fluid  approximaitons  the  model  predicted  propagation  speeds  of  the 
moving  bow  shock  to  within  five  percent  of  analytical  solutions.  Propagation  speeds  of  the  moving  magnetopause  were  also 
determined  to  within  five  percent. 


14.  SUBJECT  TERMS  - - - - 

Bow  Shock,  Magnetohydrodynamics,  Shock  Waves,  Plasma  Waves,  Space  Environment 


15.  NUMBER  OF  PAGES 

92.  ^ 

16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


18.  SECURITY  CLASSIFICATION 

OF  THIS  PAGE 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

20.  LIMITATION  OF 
ABSTRACT 

UNCLASSIFIED 

UNCLASSIFIED 

UL 

Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std,  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Dct  94 


